1
$\begingroup$

I have data consisting of soil measurements taken over 97 days at 6 time points. Soil was treated with 4 amendments (Alfalfa, Compost, Compost+Alfalfa, Control) for 288 samples (12 replications for each treatment and day).

enter image description here

My question is, what would be the model design for testing if treatment is significant for inorganic N concentrations?

I have tried lmer, I'm looking for a model with treatment and day as fixed effect and replication as a random factor.

inc.model.data <- lmer(inorg_N ~ treatment + day + (day | replication), data = mydata)
anova(inc.model.data)

mydata <- data.frame(
  treatment = c('Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Alfalfa','Alfalfa','Compost','Compost','CompAlfa','CompAlfa','Control','Control'),
  replication = c(1,2,3,4,5,6,7,8,9,10,11,12,3,4,5,6,7,8,9,10,11,12,3,4,5,6,7,8,9,10,11,12,3,4,5,6,7,8,9,10,11,12,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,1,2,1,2,1,2),
  day = c(7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,7,7,7,7,7,7,14,14),
  inorg_N = c(4.603,4.723,5.535333333,4.232333333,3.832333333,4.768333333,3.851333333,4.889333333,4.738333333,3.863333333,4.902333333,5.169333333,3.906333333,4.355333333,4.377333333,6.690333333,10.251333333,6.316333333,5.701333333,4.113333333,2.525333333,6.380333333,2.797333333,3.252333333,2.632333333,2.313333333,2.080333333,2.783333333,3.123333333,2.439333333,2.543333333,3.471333333,1.010333333,1.318333333,1.223333333,2.116333333,1.594333333,1.461333333,1.308333333,1.423333333,1.324333333,2.057333333,6.486666666,6.784666666,6.250666666,6.839666666,6.873666666,6.787666666,6.872666666,6.522666666,6.285666666,5.437666666,8.448666666,6.598666666,9.025666666,8.178666666,8.780666666,7.434666666,6.365666666,7.965666666,9.037666666,7.208666666,8.092666666,9.275666666,1.168666666,1.047666666,1.188666666,0.941666666,0.995666666,0.992666666,0.927666666,0.778666666,0.990666666,1.386666666,0.852666666,1.141666666,3.147666666,2.606666666,2.497666666,2.966666666,3.391666666,2.469666666,2.809666666,3.669666666,3.302666666,2.800666666,3.132666666,2.775666666,6.803062167,7.794062167,7.382062167,7.368062167,6.842062167,6.864062167,7.146062167,7.147062167,7.229062167,2.875062167,7.351062167,6.772062167,13.03806217,11.91006217,11.89206217,10.010062167,12.16006217,11.54106217,13.04206217,14.04306217,14.20506217,12.38106217,13.40206217,7.396062167,0.272062167,0.353062167,0.195062167,0.248062167,0.248062167,0.234062167,0.147062167,0.233062167,0.322062167,0.272062167,0.249062167,0.289062167,2.963062167,5.369062167,5.073062167,5.847062167,5.285062167,4.861062167,4.810062167,4.555062167,4.845062167,4.483062167,4.884062167,7.118062167,11.278215667,11.313215667,11.671215667,11.321215667,12.838215667,12.130215667,11.766215667,11.097215667,10.614215667,12.219215667,12.138215667,11.008215667,20.674215667,21.464215667,22.664215667,20.254215667,20.888215667,20.868215667,17.216215667,19.239215667,23.825215667,22.215215667,22.346215667,21.827215667,0.575215667,0.532215667,0.606215667,0.545215667,0.576215667,0.617215667,0.297215667,0.150215667,0.214215667,0.536215667,0.425215667,0.284215667,9.948215667,12.201215667,12.589215667,10.008215667,7.073215667,10.187215667,9.115215667,10.008215667,6.531215667,9.479215667,9.187215667,8.221215667,14.479837663,14.416837663,14.966837663,14.769837663,15.413837663,13.682837663,15.670837663,15.145837663,15.666837663,15.354837663,14.775837663,14.861837663,37.703837663,21.597837663,30.926837663,20.280837663,26.416837663,27.533837663,30.810837663,27.655837663,29.200837663,29.892837663,30.099837663,28.237837663,1.248837666,0.421837666,0.538837666,0.437837666,0.658837666,0.347837666,0.481837666,0.513837666,0.388837666,0.337837666,0.378837666,0.437837666,15.126837663,14.281837663,16.023837663,14.794837663,14.316837663,14.893837663,14.239837663,14.653837663,13.042837663,14.892837663,14.357837663,15.450837663,18.91,22.177,21.172,20.984,22.196,20.368,20.907,22.275,22.641,20.167,20.72,23.238,37.332,36.961,38.318,38.344,38.363,37.747,39.55,37.187,37.71,40.177,37.913,37.617,7.425,8.652666667,5.726,1.368,5.909666667,5.788,6.723666667,8.648666667,7.614666667,11.567,6.95,8.468,22.847666667,22.203666667,24.084666667,22.791666667,23.563666667,23.429666667,22.607666667,22.255666667,25.367666667,23.601666667,21.869666667,22.957666667,6.281333333,5.662333333,2.796333333,3.004333333,1.399333333,1.432333333,5.910333333,6.390666666))
mydata$day <- as.factor(mydata$day)

inc.model.data <- lmer(inorg_N ~ treatment + day + (day | replication), data = mydata)
anova(inc.model.data)`

enter image description here

Is this the correct approach?

An image of the experiment, there were 48 of these large jars with 6 microcosms.

Large jar with microcosms of soil amended with alfalfa

Adding summary(inc.model.data) output as per comments: enter image description here

$\endgroup$
4
  • $\begingroup$ That looks like a pretty standard mixed model for what you want to do. You are allowing the intercept to vary by the initial soil measurement within group. You then let the slope vary between plants within treatment group. Though you might want to use summary (inc.model.data) to get your beta coefficients and to see if there is any correlation between starting initial reading and day to day gain. $\endgroup$
    – JWH2006
    Jun 11 '18 at 19:21
  • $\begingroup$ Jared, can you explain your design a bit more? Do you start out with 4 batches of soil exposed to 4 treatments and then you take 12 samples per day from each batch? Can you assume those samples are independent of each other in space/time or not? $\endgroup$ Jun 12 '18 at 14:08
  • $\begingroup$ @IsabellaGhement, Gladly $\endgroup$ Jun 12 '18 at 17:40
  • 2
    $\begingroup$ The experiment was performed in an 30 C incubator. There were 336 "microcosms" of soil plus treatment. For each treatment, there were 12 large jars with 6 microcosms inside. Microcosms are the samples (I neglected to add this column to my data in the initial post), a sample was removed form the 12 large jars fore each treatment on every sampling day (this is where the 12 replicates comes in) generating 48 samples. I hope this is clear. I believe this would indicate independence...but I guess that depends on what you define independence as. Thanks again! $\endgroup$ Jun 12 '18 at 17:47
1
$\begingroup$

The fit of your model is suboptimal, unfortunately (inspect the residuals and fitted values). Here is what I would do in your situation. It is by no means the only possibility, but I think it works satisfactorily.

Here is what I have done:

  • Fitted a mixed effects model with an interaction between treatment and day for the fixed effects while dropping the random slope for day. Because the (conditional) residuals exhibited heteroscedasticity, I allowed the residual variance to vary for each combination of day and treatment.
  • Plotted the data with the predicted effects and confidence intervals as well as the original means.
  • Calculated and plotted the marginal differences between the treatments and the control for each day with 95% confidence intervals (adjusted for multiple comparisons).

Raw data with means and predicted means plus CIs

The blue line denotes the empirical means of the data. The black points with the vertical lines are the predicted means and corresponding 95% confidence intervals.

The contrasts are plotted here:

Plotted contrasts

The y-axis denotes the difference between the treatment and the control. Positive values imply that the treatment was higher compared to the control, whereas negative values imply the opposite. The vertical lines are 95% confidence intervals (adjusted for multiple comparisons).

It's clearly visible that Alfalfa exhibits higher values compared to the control for all days except day 7. CompAlfa has significantly lower values compared to the control for days 7, 14, 21 and 35. On day 49, the treatments are comparable and on day 97, CompAlfa has significantly higher values compared to the control. Compost has lower values than the control for all days.

Here is the R code I used. It contains the residual checks and the exact calculations of the contrasts.

library(nlme)
library(ggplot2)
library(emmeans)

mydata <- data.frame(
  treatment = c('Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Control','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Alfalfa','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','Compost','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','CompAlfa','Alfalfa','Alfalfa','Compost','Compost','CompAlfa','CompAlfa','Control','Control'),
  replication = c(1,2,3,4,5,6,7,8,9,10,11,12,3,4,5,6,7,8,9,10,11,12,3,4,5,6,7,8,9,10,11,12,3,4,5,6,7,8,9,10,11,12,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,3,4,5,6,7,8,9,10,11,12,1,2,1,2,1,2,1,2),
  day = c(7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,97,7,7,7,7,7,7,14,14),
  inorg_N = c(4.603,4.723,5.535333333,4.232333333,3.832333333,4.768333333,3.851333333,4.889333333,4.738333333,3.863333333,4.902333333,5.169333333,3.906333333,4.355333333,4.377333333,6.690333333,10.251333333,6.316333333,5.701333333,4.113333333,2.525333333,6.380333333,2.797333333,3.252333333,2.632333333,2.313333333,2.080333333,2.783333333,3.123333333,2.439333333,2.543333333,3.471333333,1.010333333,1.318333333,1.223333333,2.116333333,1.594333333,1.461333333,1.308333333,1.423333333,1.324333333,2.057333333,6.486666666,6.784666666,6.250666666,6.839666666,6.873666666,6.787666666,6.872666666,6.522666666,6.285666666,5.437666666,8.448666666,6.598666666,9.025666666,8.178666666,8.780666666,7.434666666,6.365666666,7.965666666,9.037666666,7.208666666,8.092666666,9.275666666,1.168666666,1.047666666,1.188666666,0.941666666,0.995666666,0.992666666,0.927666666,0.778666666,0.990666666,1.386666666,0.852666666,1.141666666,3.147666666,2.606666666,2.497666666,2.966666666,3.391666666,2.469666666,2.809666666,3.669666666,3.302666666,2.800666666,3.132666666,2.775666666,6.803062167,7.794062167,7.382062167,7.368062167,6.842062167,6.864062167,7.146062167,7.147062167,7.229062167,2.875062167,7.351062167,6.772062167,13.03806217,11.91006217,11.89206217,10.010062167,12.16006217,11.54106217,13.04206217,14.04306217,14.20506217,12.38106217,13.40206217,7.396062167,0.272062167,0.353062167,0.195062167,0.248062167,0.248062167,0.234062167,0.147062167,0.233062167,0.322062167,0.272062167,0.249062167,0.289062167,2.963062167,5.369062167,5.073062167,5.847062167,5.285062167,4.861062167,4.810062167,4.555062167,4.845062167,4.483062167,4.884062167,7.118062167,11.278215667,11.313215667,11.671215667,11.321215667,12.838215667,12.130215667,11.766215667,11.097215667,10.614215667,12.219215667,12.138215667,11.008215667,20.674215667,21.464215667,22.664215667,20.254215667,20.888215667,20.868215667,17.216215667,19.239215667,23.825215667,22.215215667,22.346215667,21.827215667,0.575215667,0.532215667,0.606215667,0.545215667,0.576215667,0.617215667,0.297215667,0.150215667,0.214215667,0.536215667,0.425215667,0.284215667,9.948215667,12.201215667,12.589215667,10.008215667,7.073215667,10.187215667,9.115215667,10.008215667,6.531215667,9.479215667,9.187215667,8.221215667,14.479837663,14.416837663,14.966837663,14.769837663,15.413837663,13.682837663,15.670837663,15.145837663,15.666837663,15.354837663,14.775837663,14.861837663,37.703837663,21.597837663,30.926837663,20.280837663,26.416837663,27.533837663,30.810837663,27.655837663,29.200837663,29.892837663,30.099837663,28.237837663,1.248837666,0.421837666,0.538837666,0.437837666,0.658837666,0.347837666,0.481837666,0.513837666,0.388837666,0.337837666,0.378837666,0.437837666,15.126837663,14.281837663,16.023837663,14.794837663,14.316837663,14.893837663,14.239837663,14.653837663,13.042837663,14.892837663,14.357837663,15.450837663,18.91,22.177,21.172,20.984,22.196,20.368,20.907,22.275,22.641,20.167,20.72,23.238,37.332,36.961,38.318,38.344,38.363,37.747,39.55,37.187,37.71,40.177,37.913,37.617,7.425,8.652666667,5.726,1.368,5.909666667,5.788,6.723666667,8.648666667,7.614666667,11.567,6.95,8.468,22.847666667,22.203666667,24.084666667,22.791666667,23.563666667,23.429666667,22.607666667,22.255666667,25.367666667,23.601666667,21.869666667,22.957666667,6.281333333,5.662333333,2.796333333,3.004333333,1.399333333,1.432333333,5.910333333,6.390666666))
mydata$day <- as.factor(mydata$day)

mydata$treatment <- relevel(mydata$treatment, ref = "Control")

inc.model.data <- lme(inorg_N~treatment * day, random=~1|replication
                      , data = mydata
                      , weights = varIdent(form= ~1|day*treatment)
                      , control = lmeControl(opt = "optim", msVerbose = TRUE))

summary(inc.model.data)
plot(inc.model.data) # Inspect residuals  
qqnorm(inc.model.data, ~resid(., type = "p")) # Q-Q-Plot of residuals

# Calculate and plot marginal means

em <- emmeans(inc.model.data, c("day", "treatment"), data = mydata)

sum_em <- summary(em)

theme_set(theme_bw())
p <- ggplot(data = mydata, aes(x = day, y = inorg_N)) +
  geom_point(aes(colour = treatment), size = 4) +
  stat_summary(aes(group = treatment), fun.y = mean,  geom = "line", size = 2, colour = "steelblue") +
  geom_pointrange(size = 1, pch = 1, data = sum_em, aes(x = day, y = emmean, ymin = lower.CL, ymax = upper.CL, group = treatment)) +
  xlab("Day") +
  ylab("Inorganic N") +
  facet_wrap(~treatment) +
  theme(
    axis.title.y=element_text(colour = "black", size = 17, hjust = 0.5, margin=margin(0,12,0,0)),
    axis.title.x=element_text(colour = "black", size = 17),
    axis.text.x=element_text(colour = "black", size=15),
    axis.text.y=element_text(colour = "black", size=15),
    legend.position="none",
    legend.text=element_text(size=12.5),
    legend.key=element_blank(),
    plot.title = element_text(face = "bold"),
    strip.text.x=element_text(size=15)
  )

p

# Calculate and plot the contrasts

em2 <- emmeans(inc.model.data, c("treatment", "day"), data = mydata)

lambdas <- list(
  "Alfalfa - Control" = c(-1, 1, rep(0, 24 - 2))
  , "CompAlfa - Control" = c(-1, 0, 1, rep(0, 24 - 3))
  , "Compost - Control" = c(-1, 0, 0, 1, rep(0, 24 - 4))
  , "Alfalfa - Control" = c(rep(0, 4), -1, 1, rep(0, 24 - 6))
  , "CompAlfa - Control" = c(rep(0, 4), -1, 0, 1, rep(0, 24 - 7))
  , "Compost - Control" = c(rep(0, 4), -1, 0, 0, 1, rep(0, 24 - 8))
  , "Alfalfa - Control" = c(rep(0, 8), -1, 1, rep(0, 24 - 10))
  , "CompAlfa - Control" = c(rep(0, 8), -1, 0, 1, rep(0, 24 - 11))
  , "Compost - Control" = c(rep(0, 8), -1, 0, 0, 1, rep(0, 24 - 12))
  , "Alfalfa - Control" = c(rep(0, 12), -1, 1, rep(0, 24 - 14))
  , "CompAlfa - Control" = c(rep(0, 12), -1, 0, 1, rep(0, 24 - 15))
  , "Compost - Control" = c(rep(0, 12), -1, 0, 0, 1, rep(0, 24 - 16))
  , "Alfalfa - Control" = c(rep(0, 16), -1, 1, rep(0, 24 - 18))
  , "CompAlfa - Control" = c(rep(0, 16), -1, 0, 1, rep(0, 24 - 19))
  , "Compost - Control" = c(rep(0, 16), -1, 0, 0, 1, rep(0, 24 - 20))
  , "Alfalfa - Control" = c(rep(0, 20), -1, 1, rep(0, 24 - 22))
  , "CompAlfa - Control" = c(rep(0, 20), -1, 0, 1, rep(0, 24 - 23))
  , "Compost - Control" = c(rep(0, 20), -1, 0, 0, 1)
)

sum_em2 <- summary(contrast(em2, lambdas), infer = c(TRUE, TRUE), adjust = "mvt")
sum_em2

sum_em2$day <- factor(rep(c(7, 14, 21, 35, 49, 97), each = 3))

theme_set(theme_bw())
p <- ggplot(data = sum_em2, aes(x = day, y = estimate)) +
  geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL, group = contrast, colour = contrast), size = 0.7) +
  geom_hline(aes(yintercept = 0), linetype = 2) +
  xlab("Day") +
  ylab("Difference Inorganic N") +
  scale_y_continuous(breaks = seq(-100, 100, 2.5)) +
  facet_wrap(~contrast) +
  theme(
    axis.title.y=element_text(colour = "black", size = 17, hjust = 0.5, margin=margin(0,12,0,0)),
    axis.title.x=element_text(colour = "black", size = 17),
    axis.text.x=element_text(colour = "black", size=15),
    axis.text.y=element_text(colour = "black", size=15),
    legend.position="none",
    legend.text=element_text(size=12.5),
    legend.key=element_blank(),
    plot.title = element_text(face = "bold"),
    strip.text.x=element_text(size=15)
  )

p
$\endgroup$
3
  • $\begingroup$ Wow, thank you! I will inspect your answer carefully, I really do appreciate your effort to help me understand this. $\endgroup$ Jun 12 '18 at 19:57
  • $\begingroup$ did you use the Q-Q plot or the scatterplot to determine heteroscedasticity? The same model without allowing the combinations of treatments and days to vary, the Q-Q plot is more S shaped instead of linear. This indicates heteroscedasticity? Thanks again. $\endgroup$ Jun 12 '18 at 21:46
  • $\begingroup$ @JaredFlater I used plot(inc.model.data) to plot the residuals against the fitted values and to check heteroscedasticity. I used the Q-Q-Plot to check the approximate normality of the residuals. $\endgroup$ Jun 13 '18 at 6:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.