Regression Analysis: Solutions to HW6

December 3, 2001

1. R&P 12.20

(a) Here are some important matrices:

> Bhat_round(B*9900,2)
Bhat_round(B*9900,2)
> Bhat
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
[1,] 85.34 56.63 62.37 53.36 48.41 51.98
[2,] 56.63 70.69 58.90 44.75 34.65 38.91
[3,] 62.37 58.90 99.20 73.06 51.28 50.09
[4,] 53.36 44.75 73.06 96.33 62.27 53.16
[5,] 48.41 34.65 51.28 62.27 75.74 58.71
[6,] 51.98 38.91 50.09 53.16 58.71 77.81
> V_matrix(0,nrow=36,ncol=36)
> V[1:6,1:6]_Bhat/4
> V[7:12,7:12]_Bhat/16
> V[13:18,13:18]_Bhat/42
> V[19:24,19:24]_Bhat/26
> V[25:30,25:30]_Bhat/9
> V[31:36,31:36]_Bhat/3
> beta_solve(t(x)%*%solve(V)%*%x)%*%t(x)%*%solve(V)%*%y
> beta
          [,1]
[1,] 0.7124256
> RSS_sum((y-beta*x)^2)
> TSS_sum(y^2)

The resulting equation is:

Ŷ_{ij} = 0.7124256X_{j}

ANOVA Table:

Source  df SS  MS
Model  1 1.582259  1.582259
Residuals 35 0.06469843 0.001848527
Total  36 1.646957

(b) Assuming common variance and common covariance in B

> mean(diag(Bhat))
[1] 84.185
> diag(Bhat)

[1] 85.34 70.69 99.20 96.33 75.74 77.81
> Bdiag_diag(c(diag(Bhat)))
> B.star_matrix(mean(Bhat-Bdiag),6,6)
> diag(B.start)_mean(diag(Bhat))
> diag(B.star)_mean(diag(Bhat))
> B.star
> B.star
         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
[1,] 84.18300 53.23533 53.23533 53.23533 53.23533 53.23533
[2,] 53.23533 84.18300 53.23533 53.23533 53.23533 53.23533
[3,] 53.23533 53.23533 84.18300 53.23533 53.23533 53.23533
[4,] 53.23533 53.23533 53.23533 84.18300 53.23533 53.23533
[5,] 53.23533 53.23533 53.23533 53.23533 84.18300 53.23533
[6,] 53.23533 53.23533 53.23533 53.23533 53.23533 84.18300
function()
{
        V.star <- matrix(0, nrow = 36, ncol = 36)
        V.star[1:6, 1:6] <- B.star/4
        V.star[7:12, 7:12] <- B.star/16
        V.star[13:18, 13:18] <- B.star/42
        V.star[19:24, 19:24] <- B.star/26
        V.star[25:30, 25:30] <- B.star/9
        V.star[31:36, 31:36] <- B.star/3
        beta.star <- solve(t(x) %*% solve(V.star) %*% x) %*% t(x) %*% solve(
                V.star) %*% y
        RSS <- sum((y - beta * x)^2)
        TSS <- sum(y^2)
        RegSS <- TSS - RSS
        MSE <- RSS/35
        return(V.star, beta.star, RSS, TSS, RegSS, MSE)
}

The resulting equation is:

Ŷ_{ij} = 0.7221587X_{j}

ANOVA Table:

Source  df SS  MS
Model  1 1.581018  1.581018
Residuals 35 0.06593967 0.001883991
Total  36 1.646957

2. Predator/Prey

y = α(1 − e^{−βx})

> summary(model.nl)
Formula: ratio ~ alpha * (1 - exp( - beta * density))
Parameters:
       Value Std. Error  t value
alpha 10.7732000  0.8509860 12.65960

Figure 1: Negative Exponential Model

 beta  0.0554244  0.0115126  4.81422
Residual standard error: 1.74606 on 24 degrees of freedom

Model Diagnostics: It is apparent that the residuals have increasing variances.

Figure 2: Model Diagnostics