Assignment 3_Heng Yan Ying
Heng Yan Ying
2026-10-08
Load packages
library(ggplot2)
library(tidyr)
library(dplyr) # mutate function
library(ggpp) # offset scatter plot on box plot
library(car) # type 2 ANOVA
library(emmeans) # stats
library(lme4) # for lmer
library(lmerTest)
library(readxl) # read excelQuestion 1: Interactions and Mixed-effects Modelling
Load and review data
lifespan <- read.csv("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 3/Assignment 3 data/lifespan.csv", header = TRUE)
# Check data structure and summary
str(lifespan)## 'data.frame': 1049 obs. of 8 variables:
## $ year : int 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
## $ Male.ID : chr "117" "99" "104" "109" ...
## $ IR.treatment: chr "20Gy" "20Gy" "Ctrl" "Ctrl" ...
## $ population : chr "C1" "C1" "C1" "C1" ...
## $ origin : chr "A30" "A30" "A30" "A30" ...
## $ date : chr "28-Dec" "28-Dec" "28-Dec" "28-Dec" ...
## $ lifespan : int 11 12 14 14 14 14 14 14 14 16 ...
## $ line : int 1 1 1 1 1 1 1 1 1 1 ...summary(lifespan)## year Male.ID IR.treatment population
## Min. :2015 Length:1049 Length:1049 Length:1049
## 1st Qu.:2015 Class :character Class :character Class :character
## Median :2016 Mode :character Mode :character Mode :character
## Mean :2016
## 3rd Qu.:2016
## Max. :2016
## origin date lifespan line
## Length:1049 Length:1049 Min. : 5.0 Min. :1.000
## Class :character Class :character 1st Qu.:12.0 1st Qu.:3.000
## Mode :character Mode :character Median :15.0 Median :5.000
## Mean :14.6 Mean :5.214
## 3rd Qu.:17.0 3rd Qu.:7.000
## Max. :29.0 Max. :9.000# Convert year to factor variable
lifespan$year <- as.factor(lifespan$year)
str(lifespan)## 'data.frame': 1049 obs. of 8 variables:
## $ year : Factor w/ 2 levels "2015","2016": 1 1 1 1 1 1 1 1 1 1 ...
## $ Male.ID : chr "117" "99" "104" "109" ...
## $ IR.treatment: chr "20Gy" "20Gy" "Ctrl" "Ctrl" ...
## $ population : chr "C1" "C1" "C1" "C1" ...
## $ origin : chr "A30" "A30" "A30" "A30" ...
## $ date : chr "28-Dec" "28-Dec" "28-Dec" "28-Dec" ...
## $ lifespan : int 11 12 14 14 14 14 14 14 14 16 ...
## $ line : int 1 1 1 1 1 1 1 1 1 1 ...Q1A
Draw out the experiment’s design and identify nested/crossed structures in the data.
Fig 1A. Experimental design diagram
Q1B
Plot the raw data to show how the lifespan of beetles varies with origin, population, and year.
lifespan %>%
ggplot(aes(x = origin, y = lifespan)) +
geom_boxplot(aes(fill = population),
width = 0.4,
outlier.shape = NA, # outliers are shown by the points
position = position_dodge(width = 0.6), # distance between boxes
alpha = 1) +
geom_point(aes(colour = population),
position = position_jitterdodge(jitter.width = 0.2,
jitter.height = 0,
dodge.width = 0.6), # match box width
alpha = 0.3, # make points more transparent
size = 0.9) +
facet_wrap(~year) + # separate both years
labs(x = "Origin",
y = "Lifespan of beetles",
fill = "Population",
colour = "Population") +
theme_minimal()
Fig 1B. Raincloud plot of lifespan of beetles against origin, separated by year and treatment.
Q1C
Plot an interaction plot between origin and year.
interaction.plot(lifespan$origin, lifespan$year, lifespan$lifespan,
xlab = "Origin",
ylab = "Lifespan of beetle",
trace.label = "Year") # label legend
Fig 1C. Interaction plot between origin and year.
Q1D
Use lmer to fit an appropriate random intercept model. (Hint: Think about what you should be including as fixed and random effects.)
# origin and year as fixed effects with * used to include interaction, population as random effect and varying intercept
lifespan_lmm <- lmer(lifespan ~ origin * year + (1 | population),
data = lifespan)Q1E
To test the significance of your fixed effects, refit the model with REML = FALSE specified as an argument, and then fit a further model without the interaction term of interest. Use anova(model1, model2, refit = FALSE) to compare the models. Do you need to keep the interaction term or not?
# Refit the model with REML = FALSE
lifespan_lmm <- lmer(lifespan ~ origin * year + (1 | population),
data = lifespan,
REML = FALSE)
# Refit the model without interaction term of interest
lifespan_lmm2 <- lmer(lifespan ~ origin + year + (1 | population),
data = lifespan,
REML = FALSE)
# Compare the models with anova
anova(lifespan_lmm, lifespan_lmm2, refit = FALSE)## Data: lifespan
## Models:
## lifespan_lmm2: lifespan ~ origin + year + (1 | population)
## lifespan_lmm: lifespan ~ origin * year + (1 | population)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## lifespan_lmm2 6 5520.2 5549.9 -2754.1 5508.2
## lifespan_lmm 8 5495.7 5535.4 -2739.8 5479.7 28.479 2 6.543e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a statistically significant interaction between origin and year (X² = 28.5, p-value < 0.001), hence the interaction term should be kept in the model.
Q1F
Refit your selected model using REML = TRUE.
# Refit the model (incl interaction) with REML = TRUE
lifespan_lmm <- lmer(lifespan ~ origin * year + (1 | population),
data = lifespan,
REML = TRUE)Q1G
Look at the summary table (i.e. summary(model)). What does it tell us about the importance of the random effect?
summary(lifespan_lmm)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: lifespan ~ origin * year + (1 | population)
## Data: lifespan
##
## REML criterion at convergence: 5480.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0126 -0.5977 -0.0044 0.6512 3.9868
##
## Random effects:
## Groups Name Variance Std.Dev.
## population (Intercept) 0.9284 0.9635
## Residual 10.7651 3.2810
## Number of obs: 1049, groups: population, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 15.24401 0.61645 5.35688 24.729 9.95e-07 ***
## originA36 -3.02709 0.36313 1040.17284 -8.336 2.41e-16 ***
## originP36 -1.15280 0.86670 5.23284 -1.330 0.239
## year2016 2.02874 0.35858 1039.10745 5.658 1.98e-08 ***
## originA36:year2016 -2.35309 0.51359 1039.22521 -4.582 5.17e-06 ***
## originP36:year2016 0.03744 0.49103 1039.44397 0.076 0.939
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) orgA36 orgP36 yr2016 oA36:2
## originA36 -0.316
## originP36 -0.711 0.225
## year2016 -0.318 0.541 0.226
## orgA36:2016 0.222 -0.703 -0.158 -0.698
## orgP36:2016 0.232 -0.395 -0.310 -0.730 0.510# Calculate percentage variance explained by population
pop_var <- 0.9284 / (0.9284 + 10.751) * 100
print(pop_var)## [1] 7.949038Only 7.95% of the variance that cannot be explained by fixed effects is explained by population. This is a small percentage which means population does not play a huge part in explaining variance that is not accounted for in fixed effects, hence not a significant random effect.
Q1H
Look at the ANOVA table. Report the appropriate F, dfs, and p-values, and interpret them.
anova(lifespan_lmm)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## origin 2865.28 1432.64 2 6.64 133.082 4.383e-06 ***
## year 406.95 406.95 1 1039.34 37.803 1.114e-09 ***
## origin:year 310.72 155.36 2 1039.40 14.432 6.575e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a statistically significant interaction between origin and year (F2,1039.4 = 14.4, p-value < 0.001). From Fig 1B, we can see that the lifespan of beetles with the A30 and P36 treatment groups had longer lifespans in the 2016 experiment compared to the 2015.
Q1I
Plot the model estimated means and intervals for your final model.
# Calculate model estimated means
lifespan_em <- emmeans(lifespan_lmm, ~ origin * year)
print(lifespan_em)## origin year emmean SE df lower.CL upper.CL
## A30 2015 15.2 0.616 5.25 13.7 16.8
## A36 2015 12.2 0.609 4.99 10.7 13.8
## P36 2015 14.1 0.609 5.01 12.5 15.7
## A30 2016 17.3 0.607 4.92 15.7 18.8
## A36 2016 11.9 0.620 5.38 10.3 13.5
## P36 2016 16.2 0.600 4.72 14.6 17.7
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95# Plot predicted means and CI
em_data <- as.data.frame(lifespan_em)
em_data %>%
ggplot(aes(x = factor(year),
y = emmean,
colour = origin)) +
geom_point(position = position_dodge(width = 0.4), size = 3) +
geom_errorbar(aes(ymin = lower.CL, # CI lower bound
ymax = upper.CL), # CI upper bound
position = position_dodge(width = 0.4),
width = 0.15) +
labs(x = "Year",
y = "Estimated mean lifespan",
colour = "Origin") +
theme_minimal()
Fig 1I. Estimated means of beetle lifespan and its 95% confidence interval for linear mixed effects model.
Question 2: Repeated Measures and Nested Analyses
Load and review data
beet <- read_excel("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 3/Assignment 3 data/Mixed models data.xlsx")
# Check data structure and summary
str(beet)## tibble [60 × 3] (S3: tbl_df/tbl/data.frame)
## $ TREATMNT: num [1:60] 1 1 1 1 1 1 1 1 1 1 ...
## $ PLANT : num [1:60] 1 1 1 1 1 1 1 1 1 1 ...
## $ DENSITY : num [1:60] 87 106.1 63.7 60.1 66.1 ...summary(beet)## TREATMNT PLANT DENSITY
## Min. :1.0 Min. :1 Min. : 12.30
## 1st Qu.:1.0 1st Qu.:1 1st Qu.: 79.78
## Median :1.5 Median :2 Median :100.00
## Mean :1.5 Mean :2 Mean : 98.03
## 3rd Qu.:2.0 3rd Qu.:3 3rd Qu.:119.05
## Max. :2.0 Max. :3 Max. :156.10Q2A
Plot the raw data to show how density varies with treatment and plant.
# Convert treatment and plant to factor variable
beet$PLANT <- as.factor(beet$PLANT)
beet$TREATMNT <- as.factor(beet$TREATMNT)
str(beet)## tibble [60 × 3] (S3: tbl_df/tbl/data.frame)
## $ TREATMNT: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ PLANT : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ DENSITY : num [1:60] 87 106.1 63.7 60.1 66.1 ...# Plot density against treatment, separated by plant
beet %>%
ggplot(aes(x = TREATMNT, y = DENSITY)) +
geom_boxplot(aes(fill = PLANT),
width = 0.4,
outlier.shape = NA, # outliers are shown by the points
position = position_dodge(width = 0.6), # distance between boxes
alpha = 1) +
geom_point(aes(colour = PLANT),
position = position_jitterdodge(jitter.width = 0.2,
jitter.height = 0,
dodge.width = 0.6), # match box width
alpha = 0.3, # make points more transparent
size = 0.9,
show.legend = FALSE) + # off legend for points
labs(x = "Soil Treatment",
y = "Bacterial Density",
fill = "Plant") +
theme_minimal()
Fig 2A. Raincloud plot of bacterial density against soil treatment, separated by plant.
Q2B
Fit the model DENSITY ~ TREATMNT, and interpret it on the assumption it is a reasonable analysis to perform.
beet_lm <- lm(DENSITY ~ TREATMNT,
data = beet)
anova(beet_lm)## Analysis of Variance Table
##
## Response: DENSITY
## Df Sum Sq Mean Sq F value Pr(>F)
## TREATMNT 1 9543 9543.2 12.504 0.0008059 ***
## Residuals 58 44266 763.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Soil treatment has a statistically significant effect on bacterial density on sugar beets’ leaf surfaces (F1,58 = 12.5, p-value = 0.000806) as p-value is less than 0.05, thus should be included in the model as a fixed effect. However, from Fig 2A, there is variation between medians of the plants hence PLANT should also be included in the model to test if it can explain the leftover variance. Additionally, as seen from the high degrees of freedom for the residuals despite only 6 plants being studied, this suggests pseudoreplication in the model.
Q2C
Recode PLANT in your dataset before fitting a model that treats plant as a random factor. Interpret this second analysis, emphasising any difference in the conclusions.
# Each plant is nested within soil treatment, create unique plant ID
beet <- beet %>%
mutate(plant_ID = interaction(TREATMNT, PLANT, drop = TRUE))
# Include plant_ID as random effect
beet_lmm <- lmer(DENSITY ~ TREATMNT + (1 | plant_ID),
data = beet)
summary(beet_lmm)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DENSITY ~ TREATMNT + (1 | plant_ID)
## Data: beet
##
## REML criterion at convergence: 544.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.57742 -0.50114 0.03571 0.68858 1.76106
##
## Random effects:
## Groups Name Variance Std.Dev.
## plant_ID (Intercept) 314.8 17.74
## Residual 546.1 23.37
## Number of obs: 60, groups: plant_ID, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 85.41 11.10 4.00 7.697 0.00153 **
## TREATMNT2 25.22 15.69 4.00 1.607 0.18329
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TREATMNT2 -0.707# Calculate variance explained by plant_ID
plant_ID_var <- 314.8/(314.8 + 546.1)*100
print(plant_ID_var)## [1] 36.56638# Checked fixed effects stats
anova(beet_lmm)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## TREATMNT 1410.6 1410.6 1 4 2.5831 0.1833In this model, soil treatment does not have a statistically significant effect on bacterial density on sugar beets’ leaf surfaces (F1,4 = 2.58, p-value = 0.183) as p-value is more than 0.05. The different plants explain 36.6% of the leftover variance after the variance explained by soil treatment. When plant is added as a random effect, soil treatment does not become a significant effect on bacterial density.
Q2D
Why can you not fit a model with a random slope and intercept model to this data?
A random slope cannot be fitted as each plant is nested within a treatment hence does not vary across treatments so a slope would not provide useful conclusions.
Q2E
For your fixed effects, report the appropriate F, dfs, and p-values from the ANOVA table. Report the importance of your random effects. Interpret your results.
From the ANOVA table, soil treatment does not have a statistically significant effect on bacterial density on sugar beets’ leaf surfaces (F1,4 = 2.58, p-value = 0.183) as p-value is more than 0.05. The different plants as a random effect explain 36.6% of the leftover variance after the variance explained by soil treatment.
Q2F
Plot the model estimated means and intervals for the model in 2C.
# Calculate model estimated means
beet_em <- emmeans(beet_lmm, ~ TREATMNT)
print(beet_em)## TREATMNT emmean SE df lower.CL upper.CL
## 1 85.4 11.1 4 54.6 116
## 2 110.6 11.1 4 79.8 141
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95# Plot predicted means and CI
em_data2 <- as.data.frame(beet_em)
em_data2 %>%
ggplot(aes(x = TREATMNT,
y = emmean)) +
geom_point(position = position_dodge(width = 0.4), size = 3) +
geom_errorbar(aes(ymin = lower.CL, # CI lower bound
ymax = upper.CL), # CI upper bound
position = position_dodge(width = 0.4),
width = 0.15) +
labs(x = "Soil Treatment",
y = "Estimated mean bacterial density") +
theme_minimal()
Fig 2F. Plot of estimated means and 95% confidence intervals of bacterial density for model DENSITY ~ TREATMNT + (1 | plant_ID).
Question 3: Split-plot Designs and Linear Mixed-Effects Models
Q3A
In this design, identify the fixed and random effects. Explain the difference between fixed and random effects.
The fixed effects are the variety of oats and level of nitrogen fertilizer, while the random effect is the blocks.
Fixed effects are the explanatory variables that we hypothesise would affect the response variable, oat yield, while random effects are grouping factors we are trying to control for, as it may influence the response variable but we are not specifically interested in it. For this design as there are more than 5 levels for the blocks, it can be entered as a random effect.
Q3B
Is this study’s design crossed or nested? Explain the difference between crossed and nested designs.
Oat variety and nitrogen fertilizer is crossed as each nitrogen fertilizer is applied to every oat variety. The oat variety is nested within the blocks. This is a nested design as each oat plant is of one variety receiving one nitrogen fertilizer treatment, nested within a block.
Crossed design is when all subjects experience all the levels of that effect, while nested design is when the factor appears only within a particular level of another factor.
Q3C
Write out your full model including any interactions and the random effect structure. (Hint: Think about the plot design.)
yield ~ nitro * Variety + (1 | Block/Variety)
Q3D
What are the number of replicates for?
Variety: 6 Nitrogen fertilizer: 18
Question 4: Power Analysis
Q4A
In the last exercise of the power analysis tutorial, you were asked to alter various aspects of the experimental design and explore different ways that the study can be performed. Select a final experimental design and describe it to us by filling in the values below:
mean_treated_greenhouses = 10.3
mean_control_greenhouses = 10
greenhouse_variation = 0.1
plant_variation = 0.1
number_of_treated_greenhouses = 5
number_of_control_greenhouses = 5
number_of_plants = 10
Power (based on 10,000 simulations) = 0.9771
Q4B
Provide justifications on why the above experimental design was selected. (Hint: Consider the number of plants and greenhouses that can be studied in a realistic scenario with limited resources.)
The above design gives a high power of 0.9771 despite the small number of greenhouses used as controls and for treatment (5 each), and a reasonable number of plants to sample per greenhouse (10 plants).
Q4C
Run a single simulation with your selected experimental design and analyse how tomato yield differs between treated and control greenhouses. Report and interpret your results.
# Set values for experiment
mean_treated_greenhouses = 10.3
mean_control_greenhouses = 10
greenhouse_variation = 0.1
plant_variation = 0.1
number_of_treated_greenhouses = 5
number_of_control_greenhouses = 5
number_of_plants = 10
# Calculate variation
variation = sqrt(greenhouse_variation^2 +
((plant_variation^2)/number_of_plants))
set.seed(0) # remember drawn values
# Give random values in a normal distribution
treated <- rnorm(number_of_treated_greenhouses,
mean = mean_treated_greenhouses,
sd = variation)
control <- rnorm(number_of_control_greenhouses,
mean = mean_control_greenhouses,
sd = variation)
# t test between treated and control greenhouse
t.test(treated, control, var.equal = TRUE)##
## Two Sample t-test
##
## data: treated and control
## t = 4.9739, df = 8, p-value = 0.001088
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.2094963 0.5716515
## sample estimates:
## mean of x mean of y
## 10.382931 9.992357There is a statistically significant difference between treated and control greenhouse plants (t8 = 4.97, p-value = 0.00109) as p-value is less than 0.05.
Q4D
Based on the analysis in 4C, plot the model estimated means and intervals.
# Save generated data into data frame
tomato <- data.frame(yield = c(treated, control),
group = factor(rep(c("Treated", "Control"),
times = c(number_of_treated_greenhouses,
number_of_control_greenhouses)),
levels = c("Control", "Treated")))
# Calculate model estimated means
tomato_lm <- lm(yield ~ group,
data = tomato)
tomato_em <- as.data.frame(emmeans(tomato_lm, ~ group))
# Plot model estimated means and CI
tomato_em %>%
ggplot(aes(x = group,
y = emmean)) +
geom_point(position = position_dodge(width = 0.4), size = 3) +
geom_errorbar(aes(ymin = lower.CL, # CI lower bound
ymax = upper.CL), # CI upper bound
position = position_dodge(width = 0.4),
width = 0.15) +
labs(x = "Greenhouse Treatment",
y = "Estimated Tomato Yield (kg)") +
theme_minimal()
Fig 4D. Plot of estimated tomato yield (kg) and 95% confidence intervals against treatment type.