Nonlinear Regression (Oddi et al. LFMC) paper

Facundo Oddi with input from Fernando Miguez

2020-04-23

Steps in fitting a nonlinear mixed model

Part I.

NONLINEAR MIXED-EFFECTS MODEL (M1)

Nonlinear modeling proccess

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

Random-Effects Mode

## 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.

## 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
## 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
##         dAIC df
## nl.me.1  0.0 20
## nl.re.2 48.5 8
##         dAIC df
## nl.me.2  0   19
## nl.me.1  2   20
##         dAIC df
## nl.me.3  0.0 18
## nl.me.2  0.4 19

Plotting the residuals

Evaluation of Fixed Effects

##            dAIC df
## nl.me.ml    0.0 21
## nl.me.ml.1  1.4 18
##            dAIC df
## nl.me.ml.1  0   18
## nl.me.ml.2  1   15
##            dAIC df
## nl.me.ml.2  0.0 15
## nl.me.ml.3 84.3 12
##            dAIC df
## nl.me.ml.2  0.0 15
## nl.me.ml.4 48.2 12

Final Model

This model assumes:

  • “A” and “w” vary with leaf type (fixed-effects).
  • “A” varies with plot (random-effects).
  • “m” and “s” are unique for all the plots and leaf types.
  • Residual variance depends on “leaf type”
  • There is no temporal dependence
  • Residuals following a normal distribution
## 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

Part II

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
# Residual checking:
par(mfrow=c(1,1))
plot(nl.fe.2) # heterocedasticity in the residuals

hist(resid(nl.fe.2), main="", xlab="Residuals") # slightly skew distribution   

qqnorm(resid(nl.fe.2))
qqline(resid(nl.fe.2))

# 3.2 ---- Final model ---------------------------------- 

M2 <- nl.fe.2
summary(M2)
## 
## 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

plot(ACF(l.me, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE) # temporal correlation is observed at lag 1

# 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

AIC(l.me, l.me.vm) # modeling the variance improves the model fit 
##         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

AIC(l.me.vm, l.me.vm.tm) # modeling the temporal correlation improves the model fit
##            df      AIC
## l.me.vm    13 2010.287
## l.me.vm.tm 15 1998.483
# Residual checking:
plot(l.me.vm.tm) 

plot(ACF(l.me.vm.tm, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE) 

hist(resid(l.me.vm.tm, type="normalized"), main="", xlab="Residuals") 

qqnorm(resid(l.me.vm.tm, type="normalized"))
qqline(resid(l.me.vm.tm, type="normalized"))

# 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
# 4.2 ---- Final model ---------------------------------- 

M3 <- l.me.vm.tm 
summary(M3)
## 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)