Assignment 3

Loading Libraries

Using the library() function to load the tidyverse package which will be used for data manipulation.

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Creating a Linear Regression Simulation

Creating the function simulate_plant_growth() to simulate plant growth data based on the soil moisture content. The simulation will generate 500 observations of soil moisture and plant growth rate The seed will be set as my student number in a hidden chunk of code.

simulate_plant_growth <- function(n = 500, 
                                  alpha = 2, 
                                  beta1 = 20, 
                                  rse = 8, 
                                  moisture_range = c(0.1, 5.0)) {
  
  soil_moisture <- runif(n, min = moisture_range[1], max = moisture_range[2])
  error <- rnorm(n, mean = 8, sd = rse)
  growth_rate <- alpha + beta1 * soil_moisture + error
  
  data.frame(
    soil_moisture = soil_moisture,
    growth_rate = growth_rate
  )
}

The runif() function generates random soil moisture values from a uniform distribution and rnorm() genrrates random error values from a normal distribution. The growth rate is calculated using the intercept, slope, soil moisture values and random error.

Running the simulation by calling upon the simulate_plant_growth() function and storing the results in an object called data_single for data manipulation and plotting. The head() function is then used to observe the first six rows of the dataframe to check the simulated data.

data_single <- simulate_plant_growth()

head(data_single)

##   soil_moisture growth_rate
## 1     2.5019810    52.02340
## 2     1.0383388    32.82435
## 3     2.8527429    61.71028
## 4     3.9070432    84.64274
## 5     0.8523106    23.48086
## 6     2.4398379    52.39973

Fit Linear Model

The lm() function is used to fit a simple linear regression model with plant growth as the response variable and soil moisture as the explanatory variable. The summary() function is then used to display the estimated coefficients, standard errors, signifiance tests and model fit statistics. The anova() function creates an Analysis Of Variance table with the sums of squares, degrees of freedom, mean squares and F-statistic which will be used in calculations and data analysis.

fit_single <- lm(growth_rate ~ soil_moisture, data = data_single)
summary(fit_single)

## 
## Call:
## lm(formula = growth_rate ~ soil_moisture, data = data_single)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.9871  -5.1989  -0.3607   4.6808  28.0428 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     9.8299     0.7817   12.58   <2e-16 ***
## soil_moisture  20.2425     0.2587   78.25   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.943 on 498 degrees of freedom
## Multiple R-squared:  0.9248, Adjusted R-squared:  0.9246 
## F-statistic:  6123 on 1 and 498 DF,  p-value: < 2.2e-16

anova(fit_single)

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

# Extract fitted values and residuals
data_single <- data_single %>%
  mutate(
    fitted = predict(fit_single),
    resid_model = growth_rate - fitted,         # Residuals from line (SSE)
    resid_mean  = growth_rate - mean(growth_rate) # Deviations from Y_bar (SST)
  )

The lm() function fitted the regression model and summary() and anova() displayed the statistical results. The %>% pipe was used to pass the data frame to the mutate() function which added three new columns. The predict() function calculated the rpedicted growth rates. the resid_model calulated the resiudals from the regression line. The resid_mean calculated the deviations from the overall mean.

Calculations

Calculating the sums of squares, mean square error, F-statistic, and R-squared value manually to verify my manual calculations using the results from the regression model.

Calculating the total sum of squares and the residual sum of squares.

SST <- sum(data_single$resid_mean^2)
SSE <- sum(data_single$resid_model^2)

The sum() function adds the squared deviations together.

Calculating the regression sum of squares

SSR <- SST - SSE

SSR was calculated by subtracting the residual sum of squares (SSE) from the total sum of squares (SST).

Calculating R-squared manually

# Calculate R-squared manually
R2_manual <- 1 - (SSE / SST)

R-squared was calculated by subtracting the proportion of unexplained variation (SSE/SST) from one

Calculating the mean square error (MSE) by dividing SSE by the residual degrees of freedom.

# Calculate MSE
MSE <- SSE / df.residual(fit_single)

Extracting the ANOVA table to obtain the regression degrees of freedom to calculate the mean square regression (MSR). Calculating the F-statistic by dividing MSR by MSE.

anova_results <- anova(fit_single)
df_regression <- anova_results["soil_moisture", "Df"]
MSR <- SSR / df_regression
F_value <- MSR / MSE

SR was calculated by dividing SSR by the regression degrees of freedom. The F-statistic was calculated by dividing MSR by MSE.

Extracting the R-squared value calculated by the lm() function to compare it to the manually calculated R-squared. The cat() function is then being used to display the calculated values for comparison with my manually calculated values.

# Compare manual R-squared with lm() output
R2_lm <- summary(fit_single)$r.squared

cat(sprintf("SST: %.2f\n", SST))

## SST: 417717.92

cat(sprintf("SSE: %.2f\n", SSE))

## SSE: 31420.84

cat(sprintf("SSR: %.2f\n", SSR))

## SSR: 386297.08

cat(sprintf("MSE: %.2f\n", MSE))

## MSE: 63.09

cat(sprintf("F-statistic: %.2f\n", F_value))

## F-statistic: 6122.56

cat(sprintf("Manual R-squared: %.5f\n", R2_manual))

## Manual R-squared: 0.92478

cat(sprintf("lm() R-squared: %.5f\n", R2_lm))

## lm() R-squared: 0.92478

The SST(), SSE, SSR, MSE, F-statistic and the R-squared all match my calculated values indicating my calculations were correct.

Creating a Plot of Plant Growth Rate Versus Soil Moisture Content

Using the ggplot function to create a scatter plot displaying the relationship between soil mositure content and plant growth rate. The fitted regression line shows the predicted relationship and the confidence band represnts uncertainty around the estimated regression line. The red lines represent the residuals.

data_single$fitted <- fitted(fit_single)
ggplot(data_single, aes(x = soil_moisture, y = growth_rate)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", color = "darkblue", se = TRUE) +
  geom_segment(aes(xend = soil_moisture, yend = fitted), alpha = 0.3, color = "red") +
  theme_classic() +
  labs(
    title = "Plant Growth Rate vs. Soil Moisture Content",
    x = "Soil Moisture Content", 
    y = "Growth Rate"
  )

## `geom_smooth()` using formula = 'y ~ x'