## Fixed-effects model for each group usign the selfStart function ("SSdlf")
nlsL <- nlsList(lfmc ~ SSdlf(time, A, w, m, s), data = lfmc.gd)
## Warning: 3 errors caught in nls(model, data = data, control = controlvals). The error messages and their frequencies are
##
## step factor 0.000488281 reduced below 'minFactor' of 0.000976562
## 1
## number of iterations exceeded maximum of 50
## 2
## A w m s
## GW plot 4 45.91299 28.66525 49.66535 -7.438673
## GW plot 5 NA NA NA NA
## GW plot 6 NA NA NA NA
## GE plot 1 101.64747 15.01433 26.93611 -11.631918
## SM plot 1 NA NA NA NA
## SS plot 1 393.28274 31.65032 22.06641 -25.829196
## GE plot 2 62.24181 12.84328 35.88690 -10.777527
## SM plot 2 292.52218 52.47574 25.63157 -15.322360
## SS plot 2 281.38633 52.34283 35.47981 -15.212427
## GE plot 3 72.58307 17.08315 27.11128 -5.651285
## SM plot 3 241.85763 16.58824 53.97104 -22.095716
## SS plot 3 236.84764 62.36245 36.17923 -10.585744
## Warning in (function (model, data = sys.frame(sys.parent()), fixed,
## random, : Iteration 1, LME step: nlminb() did not converge (code = 1). PORT
## message: function evaluation limit reached without convergence (9)
The previous model issues a warning. It is important to recognize that the unstructured variance-covariance meaning that 10 (4 + 3 + 2 + 1) parameters are tried to be estimated for the random effects. This is very likely to result in an over-parameterized model.
## Therefore,
nl.re.1 <- update(nl.re, random = pdDiag(A + w + m + s ~ 1))
## a random-effects model assuming zero correlation among the random-effects is fitted
## This converges easilty and this means the previous model (nl.re) is over parameterized.
fxf <- fixef(nl.re.1) # the intercepts of the previous model are obtained to be used in the next step.
## Mixed-effects model ("leaf type" is included as a fixed-effect on A, W, M, S)
nl.me <- update(nl.re.1, fixed = list(A + w + m + s ~ leaf.type),
start = c(A=fxf[1],0,0,0, w=fxf[2],0,0,0, m=fxf[3],0,0,0, s=fxf[4],0,0,0))
## Error in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Singularity in backsolve at level 0, block 1
## convergence problems
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: lfmc ~ SSdlf(time, A, w, m, s)
## Data: lfmc.gd
## Log-likelihood: -1070.651
## Fixed: list(A ~ 1, w ~ 1, m ~ 1, s ~ 1)
## A w m s
## 193.45254 21.45745 31.48273 -22.59708
##
## Random effects:
## Formula: list(A ~ 1, w ~ 1, m ~ 1, s ~ 1)
## Level: group
## Structure: Diagonal
## A w m s Residual
## StdDev: 119.3693 10.94687 8.924996 0.001456939 15.26567
##
## Number of Observations: 247
## Number of Groups: 12
## Therefore, the random-effects on "S" could be removed to achieve convergence
nl.re.2 <- update(nl.re.1, random = pdDiag(A + w + m ~ 1))
fxf <- fixef(nl.re.2)
## Now, the model including the fixed-effects of "leaf type" converges
nl.me.1 <- update(nl.re.2, fixed = list(A + w + m + s ~ leaf.type),
start = c(A=fxf[1],0,0,0, w=fxf[2],0,0,0, m=fxf[3],0,0,0, s=fxf[4],0,0,0))
## the covariance structure of the random effects is much smaller when the fixed-effects are modeled
nl.me.1
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: lfmc ~ SSdlf(time, A, w, m, s)
## Data: lfmc.gd
## Log-likelihood: -1034.375
## Fixed: list(A + w + m + s ~ leaf.type)
## A.(Intercept) A.leaf.typeGrass E
## 23.135414 53.294908
## A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
## 296.620596 263.313118
## w.(Intercept) w.leaf.typeGrass E
## 50.291577 -34.313924
## w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
## -26.732287 2.135665
## m.(Intercept) m.leaf.typeGrass E
## 51.533983 -21.920404
## m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus
## -20.647403 -18.170359
## s.(Intercept) s.leaf.typeGrass E
## 17.605794 -25.926031
## s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus
## -43.529450 -33.666226
##
## Random effects:
## Formula: list(A ~ 1, w ~ 1, m ~ 1)
## Level: group
## Structure: Diagonal
## A.(Intercept) w.(Intercept) m.(Intercept) Residual
## StdDev: 8.112813 0.001245815 2.468334 15.38607
##
## Number of Observations: 247
## Number of Groups: 12
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: lfmc ~ SSdlf(time, A, w, m, s)
## Data: lfmc.gd
## Log-likelihood: -1070.648
## Fixed: list(A ~ 1, w ~ 1, m ~ 1, s ~ 1)
## A w m s
## 193.46654 21.43199 31.48852 -22.60602
##
## Random effects:
## Formula: list(A ~ 1, w ~ 1, m ~ 1)
## Level: group
## Structure: Diagonal
## A w m Residual
## StdDev: 119.3792 10.94165 8.927279 15.26574
##
## Number of Observations: 247
## Number of Groups: 12
## ...suggesting parameters to vary with "leaf type"
AICtab(nl.re.2, nl.me.1) # which is supported by the AIC comparison
## dAIC df
## nl.me.1 0.0 20
## nl.re.2 48.5 8
## The small StdDev of "w" suggests to remove the random-effect of plot on this parameter
nl.me.2 <- update(nl.me.1, random = pdDiag(A + m ~ 1))
AICtab(nl.me.1, nl.me.2) # and the AIC comparison supports the simpler model
## dAIC df
## nl.me.2 0 19
## nl.me.1 2 20
## in addition, if the random-effects on "m" is also removed
nl.me.3 <- update(nl.me.2, random = pdDiag(A ~ 1))
AICtab(nl.me.2, nl.me.3) # the model fit is similar, supporting the simpler model
## dAIC df
## nl.me.3 0.0 18
## nl.me.2 0.4 19
# Heterocedasticity: variance modeling ----
nl.me.3.vm <- update(nl.me.3, weights = varIdent(form=~1|leaf.type))
AICtab(nl.me.3, nl.me.3.vm)
## dAIC df
## nl.me.3.vm 0.0 21
## nl.me.3 133.3 18
# Evaluation of the fixed-effects
nl.me.ml <- update(nl.me.3.vm, method ="ML")
## the model is fitted by Maximum likelihood
nl.me.ml.1 <- update(nl.me.ml, fixed = list(A + w + m ~ leaf.type, s ~ 1),
start = c(A=fxf[1],0,0,0, w=fxf[2],0,0,0, m=fxf[3],0,0,0, s=fxf[4]))
## the model is reduced removing the fixed-effect of "leaf type" on "s"
AICtab(nl.me.ml, nl.me.ml.1)
## dAIC df
## nl.me.ml 0.0 21
## nl.me.ml.1 1.4 18
## the AIC comparison suggests that the more complex model is not justified (fits are similar)
nl.me.ml.2 <- update(nl.me.ml, fixed = list(A + w ~ leaf.type, m + s ~ 1),
start = c(A=fxf[1],0,0,0, w=fxf[2],0,0,0, m=fxf[3], s=fxf[4]))
## the model is reduced removing the fixed-effect of "leaf type" on "m"
AICtab(nl.me.ml.1, nl.me.ml.2)
## dAIC df
## nl.me.ml.1 0 18
## nl.me.ml.2 1 15
## the more complex model is not justified
nl.me.ml.3 <- update(nl.me.ml, fixed = list(A ~ leaf.type, w + m + s ~ 1),
start = c(A=fxf[1],0,0,0, w=fxf[2], m=fxf[3], s=fxf[4]))
## the model is reduced removing the fixed-effect of "leaf type" on "W"
AICtab(nl.me.ml.2, nl.me.ml.3)
## dAIC df
## nl.me.ml.2 0.0 15
## nl.me.ml.3 84.3 12
## the AIC comparison suggests to "w" to be modeled as a function of "leaf type"
nl.me.ml.4 <- update(nl.me.ml, fixed = list(w ~ leaf.type, A + m + s ~ 1),
start = c(A=fxf[1], w=fxf[2],0,0,0, m=fxf[3], s=fxf[4]))
## the model is reduced removing the fixed-effect of "leaf type" on "A"
AICtab(nl.me.ml.2, nl.me.ml.4)
## dAIC df
## nl.me.ml.2 0.0 15
## nl.me.ml.4 48.2 12
This model assumes:
## Nonlinear mixed-effects model fit by REML
## Model: lfmc ~ SSdlf(time, A, w, m, s)
## Data: lfmc.gd
## AIC BIC logLik
## 1931.668 1983.689 -950.834
##
## Random effects:
## Formula: A ~ 1 | group
## A.(Intercept) Residual
## StdDev: 9.408521 7.162899
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | leaf.type
## Parameter estimates:
## Grass W Grass E M. spinosum S. bracteolactus
## 1.0000000 0.8851762 3.1796110 2.9647390
## Fixed effects: list(A + w ~ leaf.type, m + s ~ 1)
## Value Std.Error DF t-value p-value
## A.(Intercept) 54.25350 5.988181 226 9.060097 0e+00
## A.leaf.typeGrass E 30.92293 8.835210 226 3.499966 6e-04
## A.leaf.typeM. spinosum 223.24504 15.915558 226 14.026844 0e+00
## A.leaf.typeS. bracteolactus 240.27654 16.857355 226 14.253513 0e+00
## w.(Intercept) 29.05581 1.712960 226 16.962334 0e+00
## w.leaf.typeGrass E -20.68938 2.592469 226 -7.980570 0e+00
## w.leaf.typeM. spinosum 31.67297 7.753176 226 4.085161 1e-04
## w.leaf.typeS. bracteolactus 26.97508 8.071111 226 3.342177 1e-03
## m 30.92881 2.065196 226 14.976214 0e+00
## s -16.15406 2.365392 226 -6.829340 0e+00
## Correlation:
## A.(In) A.l.GE A..M.s A..S.b w.(In) w.l.GE
## A.leaf.typeGrass E -0.527
## A.leaf.typeM. spinosum -0.146 0.535
## A.leaf.typeS. bracteolactus -0.115 0.537 0.746
## w.(Intercept) -0.217 -0.034 -0.200 -0.216
## w.leaf.typeGrass E -0.035 -0.283 -0.382 -0.399 -0.258
## w.leaf.typeM. spinosum -0.118 -0.227 -0.547 -0.454 0.157 0.573
## w.leaf.typeS. bracteolactus -0.129 -0.241 -0.462 -0.575 0.188 0.601
## m -0.217 -0.324 -0.633 -0.666 0.092 0.145
## s -0.246 -0.354 -0.710 -0.748 0.416 0.575
## w..M.s w..S.b m
## A.leaf.typeGrass E
## A.leaf.typeM. spinosum
## A.leaf.typeS. bracteolactus
## w.(Intercept)
## w.leaf.typeGrass E
## w.leaf.typeM. spinosum
## w.leaf.typeS. bracteolactus 0.640
## m 0.161 0.172
## s 0.702 0.752 0.544
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.27999716 -0.59657157 -0.02042043 0.57273326 3.12017921
##
## Number of Observations: 247
## Number of Groups: 12
## A.(Intercept) A.leaf.typeGrass E
## 54.25350 30.92293
## A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
## 223.24504 240.27654
## w.(Intercept) w.leaf.typeGrass E
## 29.05581 -20.68938
## w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
## 31.67297 26.97508
## m s
## 30.92881 -16.15406
## A.(Intercept)
## GW plot 4 -1.915443
## GW plot 5 -1.923198
## GW plot 6 3.838641
## GE plot 1 15.652459
## SM plot 1 6.967303
## SS plot 1 9.166925
## GE plot 2 -11.066197
## SM plot 2 -5.358910
## SS plot 2 2.442781
## GE plot 3 -4.586263
## SM plot 3 -1.608393
## SS plot 3 -11.609707
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## A.(Intercept) 42.45370 54.25350 66.05331
## A.leaf.typeGrass E 13.51301 30.92293 48.33286
## A.leaf.typeM. spinosum 191.88318 223.24504 254.60691
## A.leaf.typeS. bracteolactus 207.05884 240.27654 273.49423
## w.(Intercept) 25.68039 29.05581 32.43122
## w.leaf.typeGrass E -25.79788 -20.68938 -15.58088
## w.leaf.typeM. spinosum 16.39521 31.67297 46.95073
## w.leaf.typeS. bracteolactus 11.07083 26.97508 42.87934
## m 26.85931 30.92881 34.99831
## s -20.81510 -16.15406 -11.49302
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: group
## lower est. upper
## sd(A.(Intercept)) 5.035417 9.408521 17.57953
##
## Variance function:
## lower est. upper
## Grass E 0.6772617 0.8851762 1.156919
## M. spinosum 2.4528838 3.1796110 4.121649
## S. bracteolactus 2.2864085 2.9647390 3.844316
## attr(,"label")
## [1] "Variance function:"
##
## Within-group standard error:
## lower est. upper
## 5.970704 7.162899 8.593144
The script below is identical to the portion that was distributed as the additional supporting information with the Oddi et al. publication.
# 2.2.1 - Overall predictions (Table 3) ----
Agw <- fixef(M1)[1] # "A" for grasses in the W site (GW)
Age <- fixef(M1)[1]+fixef(M1)[2] # "A" for grasses in the E site (GE)
Asm <- fixef(M1)[1]+fixef(M1)[3] # "A" for Mullinum spinosum (SM)
Ass <- fixef(M1)[1]+fixef(M1)[4] # "A" for Senecio filaginoides (SS)
Wgw <- fixef(M1)[5] # "w" for grasses in the W site (GW)
Wge <- fixef(M1)[5]+fixef(M1)[6] # "w" for grasses in the E site (GE)
Wsm <- fixef(M1)[5]+fixef(M1)[7] # "w" for Mullinum spinosum (SM)
Wss <- fixef(M1)[5]+fixef(M1)[8] # "w" for Senecio filaginoides (SS)
M <- fixef(M1)[9] # "m" for all the leaf types
S <- fixef(M1)[10] # "s" for all the leaf types
GW <- subset(lfmc, leaf.type == "Grass W")
GE <- subset(lfmc, leaf.type == "Grass E")
SM <- subset(lfmc, leaf.type == "M. spinosum")
SS <- subset(lfmc, leaf.type == "S. bracteolactus")
par(mfrow=c(2,2), cex=0.75)
# Grasses W
plot(GW$time, GW$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses W site", col="darkgreen")
curve(((Agw-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=2, lty=1, add=T, col="darkgreen")
# Grasses E
plot(GE$time, GE$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses E site", col="green")
curve(((Age-Wge)/(1 + exp((M-x)/S))+Wge), lwd=2, lty=1, add=T, col="green")
# M. Spinosum
plot(SM$time, SM$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="M. spinosum", font.main=4, col="darkorange")
curve(((Asm-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=2, lty=1, add=T, col="darkorange")
# S. filaginoides
plot(SS$time, SS$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="S. bracteolactus", font.main=4, col="darkred")
curve(((Ass-Wss)/(1 + exp((M-x)/S))+Wss), lwd=2, lty=1, add=T, col="darkred")
# 2.2.2 - Predictions at plot level (Table 3) ----
Agw.p4 <- fixef(M1)[1]+ranef(M1)[1,1] # "A" for grasses in the plot 4 (GW)
Agw.p5 <- fixef(M1)[1]+ranef(M1)[2,1] # "A" for grasses in the plot 5 (GW)
Agw.p6 <- fixef(M1)[1]+ranef(M1)[3,1] # "A" for grasses in the plot 6 (GW)
Age.p1 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[4,1] # "A" for grasses in the plot 1 (GW)
Age.p2 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[7,1] # "A" for grasses in the plot 2 (GW)
Age.p3 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[10,1] # "A" for grasses in the plot 3 (GW)
Asm.p1 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[5,1] # "A" for M. spinosum in the plot 1 (SM)
Asm.p2 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[8,1] # "A" for M. spinosum in the plot 2 (SM)
Asm.p3 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[11,1] # "A" for M. spinosum in the plot 3 (SM)
Ass.p1 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[6,1] # "A" for S. filaginoides in the plot 1 (SM)
Ass.p2 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[9,1] # "A" for S. filaginoides in the plot 2 (SM)
Ass.p3 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[12,1] # "A" for S. filaginoides in the plot 3 (SM)
par(mfrow=c(2,2), cex=0.75) # (Figure 4)
# Grasses W site
plot(GW$time, GW$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses W site", col="darkgreen")
curve(((Agw.p4-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 4
curve(((Agw.p5-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 5
curve(((Agw.p6-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 6
# Grasses E site
plot(GE$time, GE$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses E site", col="green")
curve(((Age.p1-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 1
curve(((Age.p2-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 2
curve(((Age.p3-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 3
# M. Spinosum (E site)
plot(SM$time, SM$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="M. spinosum", font.main=4, col="darkorange")
curve(((Asm.p1-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 1
curve(((Asm.p2-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 2
curve(((Asm.p3-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 3
# S. filaginoides (E site)
plot(SS$time, SS$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="S. bracteolactus", font.main=4, col="darkred")
curve(((Ass.p1-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 1
curve(((Ass.p2-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 2
curve(((Ass.p3-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 3
# 2.2.3 - Derivatives (drying speed) ----
lfmc.GW <- expression((Agw-Wgw)/(1 + exp((M-x)/S))+Wgw) # LFMC dynamics for grasses in the W site
ds.GW <- deriv(lfmc.GW, "x", function.arg = TRUE) # Drying speed for grasses in the W site
lfmc.GE <- expression((Age-Wge)/(1 + exp((M-x)/S))+Wge) # LFMC dynamics for grasses in the E site
ds.GE <- deriv(lfmc.GE, "x", function.arg = TRUE) # Drying speed for grasses in the E site
lfmc.SM <- expression((Asm-Wsm)/(1 + exp((M-x)/S))+Wsm) # LFMC dynamics for M. spinosum
ds.SM <- deriv(lfmc.SM, "x", function.arg = TRUE) # Drying speed for M. spinosum
lfmc.SS <- expression((Ass-Wss)/(1 + exp((M-x)/S))+Wss) # LFMC dynamics for S. filaginoides
ds.SS <- deriv(lfmc.SS, "x", function.arg = TRUE) # # Drying speed for S. filaginoides
xvec <- seq(0, 100, length = 1000)
y.ds.GW <- ds.GW(xvec)
y.ds.GE <- ds.GE(xvec)
y.ds.SM <- ds.SM(xvec)
y.ds.SS <- ds.SS(xvec)
par(mfrow=c(2,2), cex=0.75) # (Figure 4)
plot(xvec, y.ds.GW, xlim=c(0,80), ylim=c(0,4),
ylab="Drying speed", xlab="Time (days)", type="n", main="Grasses in the W site")
lines(xvec, -1*(attr(y.ds.GW, "grad")), lwd=2, lty=1, col="darkgreen")
plot(xvec, y.ds.GE, xlim=c(0,80), ylim=c(0,4),
ylab="Drying speed", xlab="Time (days)", type="n", main="Grasses in the E site")
lines(xvec, -1*(attr(y.ds.GE, "grad")), lwd=2, lty=1, col="green")
plot(xvec, y.ds.SM, xlim=c(0,80), ylim=c(0,4),
ylab="Drying speed", xlab="Time (days)", type="n", main="M. spinosum", font.main=4)
lines(xvec, -1*(attr(y.ds.SM, "grad")), lwd=2, lty=1, col="darkorange")
plot(xvec, y.ds.SS, xlim=c(0,80), ylim=c(0,4),
ylab="Drying speed", xlab="Time (days)", type="n", main="S. filaginoides", font.main=4)
lines(xvec, -1*(attr(y.ds.SS, "grad")), lwd=2, lty=1, col="darkred")
# 2.3 ---- Testing differences among leaf types -------------------------------------------
# 2.3.1 - Parameter "A" ----
contrast(emmeans(M1, ~leaf.type, param = "A"), "pairwise")
## contrast estimate SE df t.ratio p.value
## Grass W - Grass E -30.9 8.84 226 -3.500 0.0031
## Grass W - M. spinosum -223.2 15.92 226 -14.027 <.0001
## Grass W - S. bracteolactus -240.3 16.86 226 -14.254 <.0001
## Grass E - M. spinosum -192.3 13.45 226 -14.294 <.0001
## Grass E - S. bracteolactus -209.4 14.22 226 -14.718 <.0001
## M. spinosum - S. bracteolactus -17.0 11.71 226 -1.454 0.4671
##
## P value adjustment: tukey method for comparing a family of 4 estimates
# .- Maximum LFMC in grasses from the W site is lower than grasses from the E site (Grass W - Grass E)
# .- Maximum LFMC in grasses is lower than shrubs (Grass W - M. spinosum)
# (Grass W - S. bracteolactus)
# (Grass E - M. spinosum)
# (Grass E - S. bracteolactus)
# .- There are not difference between the two shrub species (M. spinosum - S. bracteolactus).
# 2.3.2 - Parameter "w" ----
contrast(emmeans(M1, ~leaf.type, param = "w"), "pairwise")
## contrast estimate SE df t.ratio p.value
## Grass W - Grass E 20.7 2.59 226 7.981 <.0001
## Grass W - M. spinosum -31.7 7.75 226 -4.085 0.0004
## Grass W - S. bracteolactus -27.0 8.07 226 -3.342 0.0053
## Grass E - M. spinosum -52.4 6.62 226 -7.914 <.0001
## Grass E - S. bracteolactus -47.7 6.84 226 -6.974 <.0001
## M. spinosum - S. bracteolactus 4.7 6.72 226 0.699 0.8974
##
## P value adjustment: tukey method for comparing a family of 4 estimates
# .- Minimum LFMC in grasses from the W site is higher than grasses from the E site (Grass W - Grass E)
# .- Minimum LFMC in grasses is lower than shrubs (Grass W - M. spinosum)
# (Grass W - S. bracteolactus)
# (Grass E - M. spinosum)
# (Grass E - S. bracteolactus)
# .- There are not difference between the two shrub species (M. spinosum - S. bracteolactus).
# 3) NONLINEAR FIXED-EFFECTS MODEL (M2) #####################################################################
# Response function: lfmc = (A - w) / (1 + exp((m - time)/s))) + w
# 3.1 ---- Fit of the nonlinear fixed-effects model ----------------------------------
# 3.1.1 - Base nonlinear model ----
nl.fe.0 <- nls(lfmc ~ ((A - w) / (1 + exp((m - time)/s))) + w,
start = c(A=15, m=30, s=-17, w=5),
data = lfmc) # here, time is the only predictor variable
b <- coef(nl.fe.0) # coefficients of the base model
# 3.1.2 - Leaf type fixed effect ----
nl.fe.1 <- nls(lfmc ~ ((A[leaf.type] - w[leaf.type])/(1 + exp((m - time)/s))) + w[leaf.type],
start = list(A=rep(b[1], 4), m=25, s=-10, w=rep(b[4], 4)),
data = lfmc) # and fit is made using the base model's coefficients as start values
b1 <- coef(nl.fe.1) # coefficients of the model with leaf type and time as predictor variables
# 3.1.3 - Plot fixed effect ----
nl.fe.2 <- nls(lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group],
start = list(A=rep(c(b1[1],b1[2],b1[3],b1[4]),3), m=25, s=-10, w=rep(c(b1[7],b1[8],b1[9],b1[10]),3)),
data = lfmc) # the model is fit using "b1" as the start values
AICtab(nl.fe.0, nl.fe.1, nl.fe.2) # including leaf type and plot as predictor variables imrpoves the model fit
## dAIC df
## nl.fe.2 0.0 27
## nl.fe.1 24.6 11
## nl.fe.0 684.6 5
##
## Formula: lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group]
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A1 53.683 8.602 6.240 2.20e-09 ***
## A2 54.282 8.636 6.285 1.72e-09 ***
## A3 61.906 8.881 6.971 3.59e-11 ***
## A4 111.561 13.291 8.394 5.65e-15 ***
## A5 314.131 24.840 12.646 < 2e-16 ***
## A6 333.831 26.640 12.531 < 2e-16 ***
## A7 77.465 10.997 7.044 2.34e-11 ***
## A8 298.161 24.917 11.966 < 2e-16 ***
## A9 321.433 25.765 12.475 < 2e-16 ***
## A10 87.553 11.744 7.455 2.01e-12 ***
## A11 279.437 20.829 13.416 < 2e-16 ***
## A12 292.443 24.178 12.095 < 2e-16 ***
## m 30.413 2.896 10.504 < 2e-16 ***
## s -20.695 3.606 -5.739 3.11e-08 ***
## w1 28.369 6.537 4.340 2.17e-05 ***
## w2 27.480 6.548 4.197 3.93e-05 ***
## w3 25.962 6.633 3.914 0.000121 ***
## w4 0.890 8.173 0.109 0.913388
## w5 43.524 13.567 3.208 0.001535 **
## w6 41.223 14.429 2.857 0.004686 **
## w7 6.091 7.262 0.839 0.402466
## w8 26.610 13.604 1.956 0.051725 .
## w9 39.495 14.009 2.819 0.005252 **
## w10 2.265 7.551 0.300 0.764455
## w11 66.661 11.478 5.808 2.19e-08 ***
## w12 38.222 13.031 2.933 0.003708 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.62 on 221 degrees of freedom
##
## Number of iterations to convergence: 6
## Achieved convergence tolerance: 7.269e-06
# This model assumes:
# .- "A" and "w" vary with leaf type and plot (fixed-effects).
# .- "m" and "s" are unique for all the plots and leaf types.
# .- Variance is homogeneous
# .- There is no temporal dependence
# .- Residuals following a normal distribution
# 4 - LINEAR MIXED-EFFECTS MODEL (M3) #######################################################################
# Response function: lfmc = β0 + β1xTime + β2xGE + β3xSM + β4xSS
# β5xGExTime + β6xSMxTime + β7xSSxTime
# where GE, SM, and SS are dummy variables (0 or 1) created to indicate differences between
# the base level of "leaf type", i.e., grasses from the W site [GW], and the remaining
# levels.
# 4.1 ---- Modeling process of the linear mixed-effects model ----------------------------------
# 4.1.1 - Base linear mixed-effect model ----
l.me <- lme(lfmc ~ time * leaf.type, random = ~1 | plot, data = lfmc)
# Residual checking:
plot(l.me) # heterocedasticity in the residuals
# 4.1.2 - Variance modelling ----
l.me.vm <- update(l.me, weights = varIdent(form=~ 1 | leaf.type))
# Residual checking:
plot(l.me.vm) # variances there seem to be well modeled
## df AIC
## l.me 10 2123.178
## l.me.vm 13 2010.287
# 4.1.3 - Temporal correlation modelling ----
l.me.vm.tm <- update(l.me.vm, correlation=corARMA(p=1,q=1, form=~1))
plot(ACF(l.me.vm.tm, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE) # the temporal dependence is modeled
## df AIC
## l.me.vm 13 2010.287
## l.me.vm.tm 15 1998.483
# the residuals look fine
# 4.1.4 - Evaluation of the fixed-effects ----
l.me.ml <- update(l.me.vm.tm, method ="ML") # the model is fitted by Maximum likelihood
l.me.ml.1 <- update(l.me.ml, ~.-time:leaf.type) # the model is reduced removing the interaction "time x leaf type"
AICtab(l.me.ml, l.me.ml.1) # the AIC comparison suggests the interaction term to be important
## dAIC df
## l.me.ml 0.0 15
## l.me.ml.1 51.3 12
## Linear mixed-effects model fit by REML
## Data: lfmc
## AIC BIC logLik
## 1998.483 2050.63 -984.2416
##
## Random effects:
## Formula: ~1 | plot
## (Intercept) Residual
## StdDev: 3.750488 7.842661
##
## Correlation Structure: ARMA(1,1)
## Formula: ~1 | plot
## Parameter estimate(s):
## Phi1 Theta1
## 0.9076001 -0.7567060
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | leaf.type
## Parameter estimates:
## Grass W M. spinosum S. bracteolactus Grass E
## 1.000000 2.963571 3.123740 1.070897
## Fixed effects: lfmc ~ time * leaf.type
## Value Std.Error DF t-value p-value
## (Intercept) 50.14877 3.478211 234 14.417977 0
## time -0.27535 0.048114 234 -5.722819 0
## leaf.typeGrass E 24.40614 4.908714 234 4.972002 0
## leaf.typeM. spinosum 194.85292 8.539072 234 22.818980 0
## leaf.typeS. bracteolactus 209.78398 8.838810 234 23.734415 0
## time:leaf.typeGrass E -0.56367 0.071704 234 -7.861042 0
## time:leaf.typeM. spinosum -2.02767 0.152372 234 -13.307359 0
## time:leaf.typeS. bracteolactus -2.29338 0.160937 234 -14.250193 0
## Correlation:
## (Intr) time lf.tGE lf.M.s lf.S.b tm:.GE
## time -0.557
## leaf.typeGrass E -0.709 0.395
## leaf.typeM. spinosum -0.407 0.227 0.654
## leaf.typeS. bracteolactus -0.394 0.219 0.653 0.666
## time:leaf.typeGrass E 0.374 -0.671 -0.592 -0.430 -0.429
## time:leaf.typeM. spinosum 0.176 -0.316 -0.299 -0.742 -0.409 0.546
## time:leaf.typeS. bracteolactus 0.166 -0.299 -0.318 -0.441 -0.743 0.562
## t:.M.s
## time
## leaf.typeGrass E
## leaf.typeM. spinosum
## leaf.typeS. bracteolactus
## time:leaf.typeGrass E
## time:leaf.typeM. spinosum
## time:leaf.typeS. bracteolactus 0.567
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.3449313 -0.6438200 -0.1011963 0.6158869 3.3469739
##
## Number of Observations: 247
## Number of Groups: 6
# This model assumes:
# .- β0 vary with plots as a random-effects (random intercept model)
# .- The residual variance depends on "leaf type"
# .- Temporal correlation
# .- Residuals follow a normal distribution
# 5 - CLASICAL REGRESSION (M4) ##################################################################################
M4 <- lm(lfmc ~ time * group, data = lfmc)
summary(M4)
##
## Call:
## lm(formula = lfmc ~ time * group, data = lfmc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -49.673 -7.042 0.466 7.002 66.332
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.79275 6.44309 7.573 9.60e-13 ***
## time -0.24478 0.13305 -1.840 0.06714 .
## groupGW plot 5 0.49817 9.11190 0.055 0.95645
## groupGW plot 6 6.14984 9.11190 0.675 0.50042
## groupGE plot 1 40.32675 9.49655 4.246 3.19e-05 ***
## groupSM plot 1 212.68278 9.11190 23.341 < 2e-16 ***
## groupSS plot 1 227.42468 9.11190 24.959 < 2e-16 ***
## groupGE plot 2 14.39850 9.49655 1.516 0.13089
## groupSM plot 2 195.54780 9.11190 21.461 < 2e-16 ***
## groupSS plot 2 217.12857 9.11190 23.829 < 2e-16 ***
## groupGE plot 3 21.68513 9.49655 2.283 0.02334 *
## groupSM plot 3 189.79613 9.49655 19.986 < 2e-16 ***
## groupSS plot 3 192.53490 9.49655 20.274 < 2e-16 ***
## time:groupGW plot 5 -0.01905 0.18816 -0.101 0.91944
## time:groupGW plot 6 -0.10231 0.18816 -0.544 0.58718
## time:groupGE plot 1 -0.80129 0.19357 -4.139 4.94e-05 ***
## time:groupSM plot 1 -2.36256 0.18816 -12.556 < 2e-16 ***
## time:groupSS plot 1 -2.55767 0.18816 -13.593 < 2e-16 ***
## time:groupGE plot 2 -0.43459 0.19357 -2.245 0.02574 *
## time:groupSM plot 2 -2.34722 0.18816 -12.474 < 2e-16 ***
## time:groupSS plot 2 -2.45552 0.18816 -13.050 < 2e-16 ***
## time:groupGE plot 3 -0.56657 0.19357 -2.927 0.00378 **
## time:groupSM plot 3 -1.82099 0.19357 -9.407 < 2e-16 ***
## time:groupSS plot 3 -2.16847 0.19357 -11.202 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.51 on 223 degrees of freedom
## Multiple R-squared: 0.9586, Adjusted R-squared: 0.9544
## F-statistic: 224.7 on 23 and 223 DF, p-value: < 2.2e-16
# This model assumes:
# .- The residual variance depends on "leaf type"
# .- Temporal correlation
# .- Residuals following a normal distribution
# Residual checking:
par(mfrow=c(2,2))
plot(M4) # a clear pattern in the residuals is observed (model assumptions are not met)
# 6 - NULL MODEL (M5) ######################################################################################
M5 <- lm(lfmc ~ 1, data = lfmc)
summary(M5)
##
## Call:
## lm(formula = lfmc ~ 1, data = lfmc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -88.75 -58.05 -31.45 52.66 231.06
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 95.645 4.919 19.45 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 77.3 on 246 degrees of freedom
# 7 - MODEL COMPARISON ############################################################################
# Models:
# M1: Nonlinear mixed-effects model
# M2: Nonlinear fixed-effects model
# M3: Linear mixed-effetcs model
# M4: Lnear fixed-effects model (clasical regression)
# M5: Null model
# 7.1 ---- AIC comparison (Table 1) ----------------------------------------------
AICtab(M1, M2, M3, M4, M5,
weights = T, delta = TRUE, base=T, sort = TRUE) # the nonlinear mixed-effects model is clearly the best
## AIC dAIC df weight
## M1 1931.7 0.0 15 1
## M3 1998.5 66.8 15 <0.001
## M2 2085.2 153.6 27 <0.001
## M4 2111.0 179.3 25 <0.001
## M5 2851.7 920.1 2 <0.001
# 7.2 ---- Residual analysis -------------------------------------------
par(mfrow=c(2,2), cex = 0.65) # (Figure 5)
v1 <- c(1, 2, 3, 4)
v2 <- c("GW", "GE", "SM", "SS")
# Nonlinear mixed-effects model (M1)
rM1 <- resid(M1, type = "normalized")
fM1 <- fitted(M1)
plot(fM1, rM1, xlab="Fitted LFMC (%)", ylab="Residuals", ylim=c(-3,3))
abline(a=0, b=0, lwd=2)
boxplot(rM1 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n", ylim=c(-3,3))
axis(side=1, at=v1, labels=v2, las=1)
plot(rM1 ~ lfmc$time, ylab="Residuals", xlab="Time (days)", ylim=c(-3,3))
abline(a=0,b=0, lw=2)
qqnorm(rM1, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM1)
# Nonlinear fixed-effects model (M2)
rM2 <- resid(M2)
fM2 <- fitted(M2)
plot(fM2, rM2, xlab="Fitted LFMC (%)", ylab="Residuals")
abline(a=0, b=0, lwd=2)
boxplot(rM2 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n")
axis(side=1, at=v1, labels=v2, las=1)
plot(rM2 ~ lfmc$time, ylab="Residuals", xlab="Time (days)")
abline(a=0,b=0, lw=2)
qqnorm(rM2, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM2)
# Linear mixed-effects model (M3)
rM3 <- resid(M3, type = "normalized")
fM3 <- fitted(M3)
plot(fM3, rM3, xlab="Fitted LFMC (%)", ylab="Residuals", ylim=c(-3,3))
abline(a=0, b=0, lwd=2)
boxplot(rM3 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n", ylim=c(-3,3))
axis(side=1, at=v1, labels=v2, las=1)
plot(rM3 ~ lfmc$time, ylab="Residuals", xlab="Time (days)", ylim=c(-3,3))
abline(a=0,b=0, lw=2)
qqnorm(rM3, main="", ylab="Sample quantiles", xlab="Theoretical quantiles")
qqline(rM3)
# Clasical regression (M4)
rM4 <- resid(M4)
fM4 <- fitted(M4)
plot(fM4, rM4, xlab="Fitted LFMC (%)", ylab="Residuals")
abline(a=0, b=0, lwd=2)
boxplot(rM4 ~ lfmc$leaf.type, ylab="Residuals", xlab="", xaxt="n")
axis(side=1, at=v1, labels=v2, las=1)
plot(rM4 ~ lfmc$time, ylab="Residuals", xlab="Time (days)")
abline(a=0,b=0, lw=2)
qqnorm(rM4, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM4)