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

Improve this question
$\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, 2018 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$
    – Isabella Ghement
    Jun 12, 2018 at 14:08
  • $\begingroup$ @IsabellaGhement, Gladly $\endgroup$
    – Jared Flater
    Jun 12, 2018 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$
    – Jared Flater
    Jun 12, 2018 at 17:47

1 Answer 1

Reset to default
2
$\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
Improve this answer
$\endgroup$
5
  • $\begingroup$ Wow, thank you! I will inspect your answer carefully, I really do appreciate your effort to help me understand this. $\endgroup$
    – Jared Flater
    Jun 12, 2018 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$
    – Jared Flater
    Jun 12, 2018 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$
    – COOLSerdash
    Jun 13, 2018 at 6:42
  • $\begingroup$ Hi @COOLSerdash, Tried to recreate what you said in this answer. But I had some problems interpreting my results and also how can you calculate the area under the curve from two marginal mean ? I created this question, if you could have it a look at it. Thanks. stats.stackexchange.com/questions/570145/… $\endgroup$
    – Rosa Maria
    Apr 2, 2022 at 23:44
  • $\begingroup$ Also, @COOLSerdash, how did you get the lambads for the contrast plot ? $\endgroup$
    – Rosa Maria
    Apr 2, 2022 at 23:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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