I would like to make sure this is correct even though I think it is similar to other versions on stackoverflow but not exactly the same.
Exp Design:
- Blocks - north fields and south fields
- Treatments - reference, treat_1, treat_2
- Time as months - 3, 4, 5, 6
- Response variable is nitrate - no3
The north fields have two replicates and the south fields have 1 replicate. The replicates are 2 acre fields where we measured nitrate over time in the soil as it responded to the treatments.
Packages are:
library(tidyverse)
library(car)
library(multcompView)
library(nlme)
library(emmeans)
Below is a simplified data frame.
no3.df <- structure(list(month = c(3, 3, 3, 4, 5, 5, 5, 5, 6, 3, 3, 3,
4, 5, 5, 5, 5, 6, 3, 4, 5, 5, 5, 5, 6, 3, 5, 5, 5, 5, 6, 3, 3,
3, 4, 6, 3, 3, 3, 4, 5, 5, 5, 3, 3, 4, 5, 5, 5, 5, 6, 3, 3, 3,
4, 5, 5, 5, 5, 6, 3, 3, 3, 4, 5, 5, 5, 5, 6),
block = c("north", "north", "north", "north", "north", "north", "north", "north",
"north", "north", "north", "north", "north", "north", "north",
"north", "north", "north", "south", "south", "south", "south",
"south", "south", "south", "north", "north", "north", "north",
"north", "north", "north", "north", "north", "north", "north",
"south", "south", "south", "south", "south", "south", "south",
"north", "north", "north", "north", "north", "north", "north",
"north", "north", "north", "north", "north", "north", "north",
"north", "north", "north", "south", "south", "south", "south",
"south", "south", "south", "south", "south"),
plot = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 8, 8, 8, 8, 8,
8, 8, 3, 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9, 9, 2,
2, 2, 2, 2, 2, 2, 2, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7,
7, 7, 7, 7),
treatment = c("treat_1", "treat_1", "treat_1", "treat_1",
"treat_1", "treat_1", "treat_1", "treat_1", "treat_1", "treat_1",
"treat_1", "treat_1", "treat_1", "treat_1", "treat_1", "treat_1",
"treat_1", "treat_1", "treat_1", "treat_1", "treat_1", "treat_1",
"treat_1", "treat_1", "treat_1", "treat_2", "treat_2", "treat_2",
"treat_2", "treat_2", "treat_2", "treat_2", "treat_2", "treat_2",
"treat_2", "treat_2", "treat_2", "treat_2", "treat_2", "treat_2",
"treat_2", "treat_2", "treat_2", "reference", "reference", "reference",
"reference", "reference", "reference", "reference", "reference",
"reference", "reference", "reference", "reference", "reference",
"reference", "reference", "reference", "reference", "reference",
"reference", "reference", "reference", "reference", "reference",
"reference", "reference", "reference"),
no3 = c(36.8, 20.4925, 21.03333333, 16.33, 7.723, 1.566333333, 0.533333333, 0.189, 0.31,
25.8, 16.13333333, 24.86666667, 3.979, 1.814, 0.34635, 0.244666667,
0.247333333, 0.97675, 14.305, 11.91, 12.4, 6.79, 7.26825, 8.4615,
3.43575, 22.225, 0.3243, 0.1376, 0.6244, 0.962233333, 1.36675,
8.27, 14.96, 19.62, 44.7, 9.197, 15.6, 13.85, 17.76, 14.84, 17.8,
23.06, 12.19333333, 19.06, 22.675, 27.47, 18.295, 16.5425, 18.7375,
22.25333333, 24.63125, 21.75, 23.73333333, 13.09, 20.54, 17.1,
10.58666667, 17.5565, 20.5, 25.575, 19.8, 15.76666667, 18.25333333,
15.93, 11.89, 10.791, 22.65, 22.025, 23.93333333)),
row.names = c(NA, -69L), class = c("tbl_df", "tbl", "data.frame"))
Read in the data and made factors
no3.df <- no3.df %>%
mutate(
treatment = as.factor(treatment),
plot=as.factor(plot),
month=as.factor(month))
I am using nlme to specify the covariance/variance structure. Eventually I will try this with other covariance and variance structures and look at AIC to see what is best but for now the approach I think might work best as a AR1 matrix.
lme_fitno3.block <- lme(fixed =no3 ~ treatment * month ,
random = ~1|plot/block,
method='REML',
corr = corAR1( form= ~1|plot/block),
data = no3.df)
summary(lme_fitno3.block)
Anova(lme_fitno3.block, type="III")
The model results are"
Analysis of Deviance Table (Type III tests)
Response: no3
Chisq Df Pr(>Chisq)
(Intercept) 50.8817 1 9.810e-13 ***
treatment 1.9561 2 0.376
month 3.4219 3 0.331
treatment:month 29.7859 6 4.317e-05 ***
I take from this that there is a significant interaction of treatment and month and then do followup tests.
marginal = emmeans(lme_fitno3.block,
~ treatment:month)
plot(marginal, comparisons = TRUE)
emminteraction = emmeans(lme_fitno3.block,
pairwise ~ treatment:month,
adjust="bonferroni",
alpha=0.5)
emminteraction$contrasts
multcomp::cld(marginal,
Letters = letters,
adjust="bonferroni")
I won't post the results as they are extensive.
