Assignment 3

Use my student number to reproduce the same simulated data.

1. (a)(b) Regression summary and ANOVA table

• Simulate 500 observations using the session5’s baseline settings.
• Draw moisture first, then error, to preserve the lab’s simulation order.
• Fit a simple linear regression and display the tables.

n <- 500
soil_moisture <- runif(n, min = 0.1, max = 5)
error <- rnorm(n, mean = 8, sd = 8)
growth_rate <- 2 + 20 * soil_moisture + error
plants <- data.frame(soil_moisture, growth_rate)

model <- lm(growth_rate ~ soil_moisture, data = plants)
summary(model)

## 
## Call:
## lm(formula = growth_rate ~ soil_moisture, data = plants)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.3890  -5.3931  -0.1744   5.1136  25.1254 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     8.8676     0.7471   11.87   <2e-16 ***
## soil_moisture  20.3142     0.2534   80.16   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.151 on 498 degrees of freedom
## Multiple R-squared:  0.9281, Adjusted R-squared:  0.9279 
## F-statistic:  6425 on 1 and 498 DF,  p-value: < 2.2e-16

anova(model)

## Analysis of Variance Table
## 
## Response: growth_rate
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## soil_moisture   1 426905  426905  6424.9 < 2.2e-16 ***
## Residuals     498  33090      66                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2. Equation for the observed response

• Extract the fitted intercept and slope.
• x_i is soil moisture and y_i is plant growth rate.

b0 <- unname(coef(model)[1])
b1 <- unname(coef(model)[2])

\[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i \] \[ y_i = 8.86759 + 20.31415x_i + e_i \] Here, \(e_i\) is the fitted residual.

3. Equation for the predicted response

• Store the fitted coefficients for the prediction equation.

predicted_intercept <- b0
predicted_slope <- b1

\[ \hat{y}_i = 8.86759 + 20.31415x_i \] ## 4. Equation for F - Calculate residual, total, and explained sums of squares.

SSE <- sum(residuals(model)^2)
SST <- sum((plants$growth_rate - mean(plants$growth_rate))^2)
SSR <- SST - SSE

df_error <- df.residual(model)
F_value <- (SSR / 1) / (SSE / df_error)

\[ F = \frac{MSR}{MSE} = \frac{SSR/1}{SSE/(n-2)} \]

\[ F = \frac{426905.48073/1} {33089.74275/498} = 6424.91938 \]

5. Equation for MSE

• Divide residual variation by the residual degrees of freedom.

MSE <- SSE / df_error

\[ MSE = \frac{SSE}{n-2} \]

\[ MSE = \frac{33089.74275}{498} = 66.44527 \]

6. Equation for SST

• Calculate total variation around the sample mean.

SST <- sum((plants$growth_rate - mean(plants$growth_rate))^2)

\[ SST = \sum_{i=1}^{n}(y_i-\bar{y})^2 = SSR + SSE \]

\[ SST = 426905.48073 + 33089.74275 = 459995.22348 \]

7. Equation for r squared

• Calculate the proportion of variation explained by the model.

R_squared <- 1 - SSE / SST

\[ r^2 = SSR / SST = 1 - SSE / SST \]

\[ r^2 = 1 - 33089.74275 / 459995.22348 = 0.92807 \]

8. Does soil moisture content affect plant growth rate?

• Extract the model statistics and slope confidence interval.

model_summary <- summary(model)
f_stat <- model_summary$fstatistic
slope <- unname(coef(model)["soil_moisture"])
slope_ci <- confint(model, "soil_moisture")
p_value <- pf(f_stat[1], f_stat[2], f_stat[3], lower.tail = FALSE)

Soil moisture was associated with plant growth rate in the simulated dataset (linear regression: F(1, 498) = 6424.92, p < 0.001, R² = 0.9281, n = 500). Each one-unit increase in soil moisture was associated with an estimated 20.31-unit increase in growth rate (95% CI: 19.82 to 20.81). The model explained 92.81% of the variation in growth rate.

9. Regression plot with confidence band

• Predict mean growth and its 95% confidence interval along the moisture range.
• Set up the axes, then add the band, observations, and fitted line.

new_moisture <- data.frame(
  soil_moisture = seq(min(soil_moisture), max(soil_moisture), length.out = 200)
)
predictions <- predict(model, newdata = new_moisture,
                       interval = "confidence", level = 0.95)


plot(growth_rate ~ soil_moisture, data = plants, type = "n",
     ylim = range(growth_rate, predictions),
     xlab = "Soil moisture content", ylab = "Plant growth rate")
x <- new_moisture$soil_moisture
polygon(c(x, rev(x)),
        c(predictions[, "lwr"], rev(predictions[, "upr"])),
        col = "lightblue", border = NA)
points(plants$soil_moisture, plants$growth_rate,
       pch = 16, cex = 0.6, col = "black")
lines(x, predictions[, "fit"], col = "blue", lwd = 2)