Question 1: Interactions and Mixed-effects Modelling

The dataset lifespan.csv contains data from Berger, D., Stångberg, J., Grieshop, K., Martinossi-Allibert, I. & Arnqvist, G. (2017). Temperature effects on life-history trade-offs, germline maintenance and mutation rate under simulated climate warming. Proceedings. Biological sciences / The Royal Society, 284.

The authors were interested in the effects of warming, and of acclimatisation to warming on a variety of aspects of the biology of the bruchid beetle Callosobruchus maculatus. They exposed three lines (“preadapted lines”) to warming, with temperatures increasing over time until they reached 36ºC, and maintained another three control lines (“ambient”) at 30ª, and three other control lines at 36ºC. Individuals from all lines were then monitored for lifespan giving three treatments at either 30ºC or 36ª (controls), or 36ª (preadapted lines).

lifespan - Our response variable: how long an animal lived in days. origin - A factor which tells us whether the animal was from a control (ambient) line at 30ª (A30), a control line at 36ª (A36) or a preadapted line at 36ª (P36). population - The specific line that the individual came from year - The year that the experiment was carried out: either 2015 or 2016.

The question we would like to answer is how longevity is affected by preadaptation and temperature, and whether these effects vary with year. The 2016 data are from lines where the preadapted beetles had been exposed to the higher temperature for 15 generations longer than in 2015, so we might expect a difference. The data are grouped by population, so we need to include this as a random effect in our model.

We’ll use lmer from the lme4 package to fit our model.

1A. Draw out the experiment’s design and identify nested/crossed structures in the data.

1A. Experiment Design Structure

There is a crossed structure composing of the year and the type of exposed line while there is a nested structure between the type of exposed line and how many lines there are present for each, in this case 3. There is also further nesting structure between the number of beetle individuals lifespan monitored for each line.

1B. Plot the raw data to show how the lifespan of beetles varies with origin, population, and year.

# extract data from csv

beetles<-read.csv('C:/Users/mishr/Desktop/Y4S1/ES3307/Assignment 3/lifespan.csv')

summary(beetles)

##       year           Male.ID        IR.treatment      population  
##  Min.   :2015   Length   :1049   Length   :1049   Length   :1049  
##  1st Qu.:2015   N.unique : 804   N.unique :   2   N.unique :   6  
##  Median :2016   N.blank  :   0   N.blank  :   0   N.blank  :   0  
##  Mean   :2016   Min.nchar:   1   Min.nchar:   4   Min.nchar:   2  
##  3rd Qu.:2016   Max.nchar:   6   Max.nchar:   4   Max.nchar:   2  
##  Max.   :2016                                                     
##        origin            date         lifespan         line      
##  Length   :1049   Length   :1049   Min.   : 5.0   Min.   :1.000  
##  N.unique :   3   N.unique :  15   1st Qu.:12.0   1st Qu.:3.000  
##  N.blank  :   0   N.blank  :   0   Median :15.0   Median :5.000  
##  Min.nchar:   3   Min.nchar:   5   Mean   :14.6   Mean   :5.214  
##  Max.nchar:   3   Max.nchar:   6   3rd Qu.:17.0   3rd Qu.:7.000  
##                                    Max.   :29.0   Max.   :9.000

# plot all 3 variables with lifespan of beetles on a combined figure #
ggplot(beetles, aes(x = origin, y = lifespan, colour = population)) +
  geom_boxplot(aes(group = origin),
               colour = "grey40", alpha = 0.25,
               outlier.shape = NA) +
  geom_jitter(width = 0.15, height = 0,
              alpha = 0.65, size = 1.8) +
  facet_wrap(~ year) +
  labs(
    title = "Variation of Beetle lifespan by Origin, Population and Year",
    x = "Origin",
    y = "Lifespan (days)",
    colour = "Population"
  ) +
  theme_classic()

The beetle lifespan distributions vary across origins and the two years, with A30 (ambient line) generally showing higher median lifespan in 2016 than 2015, while A36 remains relatively low in both years. The distribution spread is the largest for A30 in both years, implying greater variation in lifespan, compared to A36 (control line), which has a narrower distribution. P36 (preadapted line) shows a similar median beetle lifespan across years with moderate variation, suggesting that lifespan patterns does differ by origin and year.

1C. Plot an interaction plot between origin and year.

ggplot(beetles, aes(x = origin, y = lifespan,
                     colour = factor(year),
                     group = factor(year))) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", linewidth = 0.8) +
  stat_summary(fun.data = mean_se, geom = "errorbar",
               width = 0.12) +
  labs(
    title = "Interaction between Origin and Year",
    x = "Origin",
    y = "Mean lifespan (days)",
    colour = "Year"
  ) +
  theme_classic()

1D. Use lmer to fit an appropriate random intercept model. (Hint: Think about what you should be including as fixed and random effects.)

model_beetles <- lmer(
  lifespan ~ origin * year + (1 | population),
  data = beetles,
)

summary(model_beetles)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: lifespan ~ origin * year + (1 | population)
##    Data: beetles
## 
## 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)    -4.073e+03  7.227e+02  1.038e+03  -5.635 2.25e-08 ***
## originA36       4.738e+03  1.035e+03  1.039e+03   4.578 5.27e-06 ***
## originP36      -7.660e+01  9.897e+02  1.039e+03  -0.077    0.938    
## year            2.029e+00  3.586e-01  1.038e+03   5.658 1.98e-08 ***
## originA36:year -2.353e+00  5.136e-01  1.039e+03  -4.582 5.17e-06 ***
## originP36:year  3.744e-02  4.910e-01  1.039e+03   0.076    0.939    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) orgA36 orgP36 year   orA36:
## originA36   -0.698                            
## originP36   -0.730  0.510                     
## year        -1.000  0.698  0.730              
## orignA36:yr  0.698 -1.000 -0.510 -0.698       
## orignP36:yr  0.730 -0.510 -1.000 -0.730  0.510
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling. 
## You may also use (g)lmerControl(autoscale = TRUE) to improve numerical stability.

The origin and year are fixed effects while population is random effect.

1E. 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?

model_int <- lmer(
  lifespan ~ origin * year + (1 | population), # interaction indicated by (*)
  data = beetles,
  REML = FALSE
)

model_no_int <- lmer(
  lifespan ~ origin + year + (1 | population),
  data = beetles,
  REML = FALSE
)

anova(model_no_int, model_int, refit = FALSE)

## Data: beetles
## Models:
## model_no_int: lifespan ~ origin + year + (1 | population)
## model_int: lifespan ~ origin * year + (1 | population)
##              npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
## model_no_int    6 5520.2 5549.9 -2754.1    5508.2                         
## model_int       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 ' ' 1

compare_models_1<-anova(model_no_int, model_int, refit = FALSE)


# Create kable table
knitr::kable(
  compare_models_1,
  digits = 3,
  caption = "Anova Results Comparing the Two Models"
)

Anova Results Comparing the Two Models
	npar	AIC	BIC	logLik	-2*log(L)	Chisq	Df	Pr(>Chisq)
model_no_int	6	5520.188	5549.921	-2754.094	5508.188	NA	NA	NA
model_int	8	5495.708	5535.353	-2739.854	5479.708	28.479	2	0

The interaction between origin and year was statistically significant (X^2 (2) =28.48, p<0.001). This implies that there is very strong statistical evidence to support that the effect of origin on beetle lifespan indeed differs between 2015 and 2016. Thus, the interaction term should be retained.

1F. Refit your selected model using REML = TRUE.

model_1f <- lmer(
  lifespan ~ origin * year + (1 | population),
  data = beetles,
  REML = TRUE
)

1G. Look at the summary table (i.e. summary(model)). What does it tell us about the importance of the random effect?

summary(model_1f)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: lifespan ~ origin * year + (1 | population)
##    Data: beetles
## 
## 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)    -4.073e+03  7.227e+02  1.038e+03  -5.635 2.25e-08 ***
## originA36       4.738e+03  1.035e+03  1.039e+03   4.578 5.27e-06 ***
## originP36      -7.660e+01  9.897e+02  1.039e+03  -0.077    0.938    
## year            2.029e+00  3.586e-01  1.038e+03   5.658 1.98e-08 ***
## originA36:year -2.353e+00  5.136e-01  1.039e+03  -4.582 5.17e-06 ***
## originP36:year  3.744e-02  4.910e-01  1.039e+03   0.076    0.939    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) orgA36 orgP36 year   orA36:
## originA36   -0.698                            
## originP36   -0.730  0.510                     
## year        -1.000  0.698  0.730              
## orignA36:yr  0.698 -1.000 -0.510 -0.698       
## orignP36:yr  0.730 -0.510 -1.000 -0.730  0.510
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling. 
## You may also use (g)lmerControl(autoscale = TRUE) to improve numerical stability.

The random effect of population had a variance of 0.93 (2 d.p.), compared with the residual variance of 10.77 (2 d.p.). This implies that variation between populations was relatively small compared with variation between individual beetles. The population random effect accounts for about 7.9% (0.9284 / (0.9284 + 10.7651)) of the total estimated variance, indicating a moderate population-level effect.

1H. Look at the ANOVA table. Report the appropriate F, dfs, and p-values, and interpret them.

anova(model_1f)

## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## origin      310.26  155.13     2 1038.9  14.410 6.713e-07 ***
## year        406.95  406.95     1 1038.9  37.803 1.114e-09 ***
## origin:year 310.72  155.36     2 1038.9  14.432 6.575e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova_models_1f<-anova(model_1f)

# Create kable table
knitr::kable(
  anova_models_1f,
  digits = 3,
  caption = "Anova Results"
)

Anova Results
	Sum Sq	Mean Sq	NumDF	DenDF	F value
origin	310.256	155.128	2	1038.884	14.410
year	406.953	406.953	1	1038.909	37.803
origin:year	310.715	155.358	2	1038.883	14.432

The ANOVA results show that there is very strong statistical evidence that effect origin on beetle lifespan is significant (F_2,1038.88 =14.41, p<0.001). Similarly, there was also a significant effect of year (F~1,1038.90~ = 37.80, p < 0.001). In a similar vein, the interaction between origin and year was also statistically significant (F_2,1038.88 = 14.43, p <0.001), implying that the effect of origin on beetle lifespan indeed differed between 2015 and 2016.

1I. Plot the model estimated means and intervals for your final model.

emm <- emmeans(model_1f, ~ origin * year)

emm_1i <- as.data.frame(emm)

ggplot(emm_1i,
       aes(x = origin, y = emmean,
           colour = factor(year),
           group = factor(year))) +
  geom_point(position = position_dodge(width = 0.15),
             size = 3) +
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.15),
    width = 0.1
  ) +
  geom_line(linewidth = 0.8) +
  labs(
    title = "Model Estimated Means and Intervals",
    x = "Origin",
    y = "Estimated mean lifespan (days)",
    colour = "Year"
  ) +
  theme_classic()

Question 2: Repeated Measures and Nested Analyses

An experiment was performed to determine the effect of two soil treatments on the density of bacteria found on the leaf surfaces of sugar beet. Six plants were grown in an environmental chamber, three in Soil type 1 and three in Soil type 2. Then samples were taken from each of the six plants. For each sample, a fixed biomass of leaf was homogenised, and a fixed volume of a designated dilution spread on a plate. The bacterial density was then assessed by counting the number of colonies for each plate.

The data were stored in variable TREATMNT (1 or 2) and PLANT is coded 1 to 3 in each treatment. The number of colonies are recorded in DENSITY. Two sets of analyses were conducted on these data (dataset Mixed models data.xlsx), to ask if the soil treatment influenced the density of bacteria on the leaf surface.

2A.Plot the raw data to show how density varies with treatment and plant.

# load .xlsx dataset
soil<-read_excel("C:/Users/mishr/Desktop/Y4S1/ES3307/Assignment 3/Mixed models data.xlsx")

# check data types
glimpse(soil)

## Rows: 60
## Columns: 3
## $ TREATMNT <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ PLANT    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3…
## $ DENSITY  <dbl> 87.0, 106.1, 63.7, 60.1, 66.1, 37.9, 101.4, 12.3, 74.7, 54.7,…

# currently all saved as dbl, meaning floating numbers/numeric values
# as treatment and plant are not numeric but rather categorical we need to use factor()

# so convert treatment and plant to factors 
soil$TREATMNT<- as.factor(soil$TREATMNT)
soil$PLANT<- as.factor(soil$PLANT)

# check 
str(soil)

## 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 boxplot for density versus plant # 
soil %>% 
  ggplot(aes(x = PLANT, y = DENSITY)) +
  geom_boxplot(fill = "purple", outlier.shape = NA, width = 0.15) + 
  geom_jitter(color = "steelblue", width = 0.2, alpha = 0.6, size = 2) +
  labs(title = "Relationship between Plant and Bacterial Density ", 
       x = "Plant", 
       y = "Number of Colonies") +
  theme_light()

# plot boxplot for density versus treatment # 
soil %>% 
  ggplot(aes(x = TREATMNT, y = DENSITY)) +
  geom_boxplot(fill = "yellow", outlier.shape = NA, width = 0.15) + 
  geom_jitter(color = "red", width = 0.2, alpha = 0.6, size = 2) +
  labs(title = "Effect of Soil Treatment on Bacterial Density ", 
       x = "Treatment", 
       y = "Number of Colonies") +
  theme_light()

From the bacterial density versus plant boxplot, we can observe that the distribution of number of colonies (bacterial density) is similar across the 3 plants, indicating likely equal variance. No particular significant outlier present. We can reasonably infer that the data is normally distributed.

From the bacterial density versus soil treatment boxplot, we can observe that Treatment 1 has a more spread out but lower distribution of points than treatment 2. Hard to tell if there is a statistical significant difference between the mean bacterial density between the two treatments just visually, better to perform further tests.

2B. Fit the model DENSITY ~ TREATMNT, and interpret it on the assumption it is a reasonable analysis to perform.

soil_m1 <- lm(DENSITY ~ TREATMNT, data = soil)
anova(soil_m1)

## 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 ' ' 1

par(mfrow = c(2, 2))
plot(soil_m1)

anova_models_2b<-anova(soil_m1)

# Create kable table
knitr::kable(
  anova_models_2b,
  digits = 3,
  caption = "Anova Results"
)

Anova Results
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
TREATMNT	1	9543.248	9543.248	12.504	0.001
Residuals	58	44265.884	763.205	NA	NA

From the model results, we can observe that there is very strong statistical evidence to show that the difference in mean bacterial density between the two soil treatments is significant ((F_1,58 = 12.50, p < 0.001)). Thus, we reject the null hypothesis as soil treatment has a statistically significant effect on bacterial density.

From the diagnostic plots, - Residuals vs Fitted: Residuals form two distinct vertical groups because there are two treatment levels, each with its own fitted mean.

Normal Q-Q: Residuals follow the straight line well, indicating that they are normally distributed.
Scale-Location: The horizontal red line indicates that homoscedasticity can be inferred.
Constant Leverage: Residuals vs Factor Levels: Both groups center near zero, with no particular influential outlier present.

2C. 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.

soil_m2 <- lmer(DENSITY ~ TREATMNT + (1 | PLANT), data = soil)
anova(soil_m2)

## Type III Analysis of Variance Table with Satterthwaite's method
##          Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## TREATMNT 9543.2  9543.2     1    56  17.816 9.006e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

soil_m2_anova <- anova(soil_m2)

# Create kable table
knitr::kable(
  soil_m2_anova,
  digits = 3,
  caption = "Anova results of Treating Plant as a Random Factor"
)

Anova results of Treating Plant as a Random Factor
	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
TREATMNT	9543.248	9543.248	1	56	17.816	0

By treating plant as a random factor, the model gives a result of F~1, 56~ = 17.82, p < 0.001, so the effect of soil treatment on bacterial density is likely still highly significant. This conclusion matches 2B, where we reject the null hypothesis which states that there is no statistical significant difference on bacterial density between the two treatments; but in this case the evidence is even stronger. F increases from 12.50 to 17.82 (2 d.p.) and the p-value decreases to about 9 × 10⁻⁵. This is because the random plant variable absorbs the variation between plants, which decreases the error term used to test treatment. The denominator df (56) is Satterthwaite-approximated, so it differs slightly from the 58 in 2B, where plant was treated as fixed effect.

2D. Why can you not fit a model with a random slope and intercept model to this data?

In order to fit a model with a random slope or intercept model, there needs to be sufficient levels, but in our case there are only 3 plants, which means too few grouping levels. Furthermore, a 2×2 covariance matrix (intercept variance, slope variance, and their correlation) is required, and 3 levels is not enough, so the fit cannot be properly accounted for.

2E. 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.

summary(soil_m2)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DENSITY ~ TREATMNT + (1 | PLANT)
##    Data: soil
## 
## REML criterion at convergence: 541
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.74839 -0.48261 -0.03142  0.62083  1.81966 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PLANT    (Intercept) 330.0    18.16   
##  Residual             535.6    23.14   
## Number of obs: 60, groups:  PLANT, 3
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   85.413     11.307  2.311   7.554   0.0112 *  
## TREATMNT2     25.223      5.976 56.000   4.221 9.01e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## TREATMNT2 -0.264

From 2D, indeed, soil treatment has a significant effect on bacterial density, F_1,56 = 17.82, p < 0.001, with treatment 2 having a higher mean of about 25.22 (2 d.p.) more colonies than treatment 1. The random effect tells us how much variance we find among levels of our grouping factor and the residual factor. The random effect, in this case the plant is important, as plant variance (330.0, SD = 18.2) is substantial relative to residual variance (535.6, SD = 23.14). Thus, by accounting for plant-to-plant variation, this reduces the error and increases the F value, compared with 2b (17.82 vs. 12.50) (2 d.p.), but the conclusion still stands.

2E. Plot the model estimated means and intervals for the model in 2C.

# Estimated marginal means for each treatment
em_m2 <- emmeans(soil_m2, ~ TREATMNT)
em_m2_df <- as.data.frame(em_m2)
em_m2_df

##  TREATMNT    emmean       SE   df lower.CL upper.CL
##  1         85.41333 11.30669 2.31 42.53150 128.2952
##  2        110.63667 11.30669 2.31 67.75483 153.5185
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

# Plot
ggplot(em_m2_df, aes(x = TREATMNT, y = emmean)) +
  geom_point(size = 3, colour = "darkred") +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                width = 0.1, colour = "darkred") +
  labs(title = "Model Estimated Means and Intervals",
       x = "Soil Treatment",
       y = "Estimated mean bacterial density (colonies)") +
  theme_light()

Question 3: Split-plot Designs and Linear Mixed-Effects Models

Above is a schematic of the experimental design for the OATS study by Yates (1935). The study site was divided into six blocks. Within each block, there are three transects, each planted with one of the three varieties of oats. Within each variety transect, there are four subplots, each with a different level of nitrogen fertilizer. The response variable is yield at the end of the growing season.

The data variables are as follows:

Block: An ordered factor with levels VI < V < III < IV < II < I
Variety: A factor with levels Golden Rain, Marvellous, Victory
nitro: A numeric vector of 4 nitrogen fertilizer concentrations
yield: A numeric vector of the yield of oats

Experimental Design, OATS study by Yates (1935)

3A. In this design, identify the fixed and random effects. Explain the difference between fixed and random effects.

Qn 3. Experimental Design Structure

In this design, the fixed effects are the variety of oats, where there are 3 different kinds, namely Golden Rain, Marvellous and Victory, the nitrogen fertilizer concentration and the variety and nitrogen fertilizer concentration interaction. Whereas, the random effects are the blocks, transect arrangement in each block and the arrangement of the different nitrogen fertilizer concentrations in each transect.

A fixed effect refers to variables that we are interested in and how they might affect our experiment or other treatments. They are not randomly assigned, unlike random effects which refer to variables that may affect the experiment but we are not particularly interested in their effects.

3B. Is this study’s design crossed or nested? Explain the difference between crossed and nested designs.

(attach image)

This study’s design is split plot where there is both crossing and nesting present. For a crossed design, all levels of one factor, in this case oat variety, is tested with all levels of nitrogen fertilizer concentration, and is repeated across all blocks. On the other hand, a nested design would mean that one factor (nitrogen fertilizer concentration, A) is nested within another factor (variety, B), where the levels of that levels of that factor (A) are only observed at specific levels of another factor (B). Therefore, variety and nitrogen fertilizer concentration are crossed treatment factors, while the transect structure is nested within blocks.

3C. Write out your full model including any interactions and the random effect structure. (Hint: Think about the plot design.)

oat_model <- lmer(yield ~ variety*nitrogen_fertilizer_con + (1|block/variety), data = oats)

In R, (*), represents interactions.

3D.What are the number of replicates for?

Variety: 6
Nitrogen fertilizer: 18 (6 x 3)

Question 4: Power Analysis

4A. 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.153
plant_variation = 0.5
number_of_treated_greenhouses = 6
number_of_control_greenhouses = 6
number_of_plants = 80
Power (based on 10,000 simulations) = 0.82 (2 d.p.)

4B. 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.)

For an experiment that would actually be carried out, one needs to critically consider what is the most apt way to optimise limited resources while still ensuring a balanced and well thought out experimental design. To ensure a balanced design, we keep the number of treated and control greenhouses the same, in this case 6 each so in total 12. The reason I increased the number of greenhouses is so as to increase the number of independent experimental replicates, which can increase statistical power. In a realistic scenario, 80 plants may be “easier” to monitor for yield changes (response variable), although the overall 80 x 12 = 960 is still relatively high, it is more manageable than 100 x 12 = 1200 plants to observe for yield. The statistical power is about 0.82 (2 d.p.), implying that the experiment has about 82% probability of detecting a statistically significant treatment effect. Hence, in a study with limited resources, these considerations may make the execution of this experiment more realistic.

4C. 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.

# Alter the values below
mean_treated_greenhouses = 10.3
mean_control_greenhouses = 10
greenhouse_variation = 0.153
plant_variation = 0.5

number_of_treated_greenhouses = 6
number_of_control_greenhouses = 6
number_of_plants = 80

# for reproducability
set.seed(123)

# Do not edit anything below here
results=vector()
correct=0
runs = 10000
variation = sqrt(greenhouse_variation^2 + ((plant_variation^2)/number_of_plants))
variation

## [1] 0.1628926

for(i in 1:runs){
  treated<-rnorm(number_of_treated_greenhouses, mean = mean_treated_greenhouses, sd = variation)
  
  control<-rnorm(number_of_control_greenhouses, mean = mean_control_greenhouses, sd = variation)
  
  
  model1<-t.test(treated, control, var.equal = TRUE)
  
  p<-model1$p.value
  if(p<0.05) correct = correct +1
  results[i]<-p
}

power = correct/runs
power

## [1] 0.8165

data.frame(treated,control)

##    treated   control
## 1 10.50285  9.934727
## 2 10.60271 10.015902
## 3 10.26553 10.227196
## 4 10.30538 10.029370
## 5 10.41466  9.995571
## 6 10.47058  9.895932

model1 <- t.test(treated, control, var.equal = TRUE)
model1

## 
##  Two Sample t-test
## 
## data:  treated and control
## t = 5.8891, df = 10, p-value = 0.0001533
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.2551885 0.5658132
## sample estimates:
## mean of x mean of y 
##  10.42695  10.01645

model1_results <- data.frame(
  Statistic = model1$statistic,
  df = model1$parameter,
  p_value = model1$p.value,
  Mean_Treated = model1$estimate[1],
  Mean_Control = model1$estimate[2]
)

# Generate the beautiful table
knitr::kable(
  model1_results,
  digits = 3,
  caption = "Results of A Single Simulation"
)

Results of A Single Simulation
	Statistic	df	p_value	Mean_Treated	Mean_Control
t	5.889	10	0	10.427	10.016

The two-sample t-test showed a statistically significant difference between the yield of the plants from the treated and controlled greenhouse groups (t(10)~ = 5.89, p < 0.001). Thus, we reject the null hypothesis that there is no true difference between the controlled and treatment groups.

4D. Based on the analysis in 4C, plot the model estimated means and intervals.

plot_4d <- data.frame(Yield = c(treated, control),
                        Group = factor(rep(c("Treated","Control"), each = 6),
                                       levels = c("Control","Treated")))
fit <- lm(Yield ~ Group, data = plot_4d)
em  <- as.data.frame(emmeans(fit, ~ Group))

ggplot(em, aes(Group, emmean)) +
  geom_jitter(data = plot_4d, aes(y = Yield), width = 0.08, alpha = 0.5) +
  geom_point(size = 4, colour = "darkgreen") +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
  labs(y = "Mean yield from each greenhouse group (kg)", x = NULL,
       title = "Model Estimated Means and Intervals") +
  theme_minimal()

ES3307 Assignment 3

Misha Rashanah

2026-10-16