Nonlinear Regression (Oddi et al. LFMC) paper

Facundo Oddi with input from Fernando Miguez

2020-06-12

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