# Steps in fitting a nonlinear mixed model

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

### 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 w..M.s
## 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  0.640
## m                           -0.217 -0.324 -0.633 -0.666  0.092  0.145  0.161
## s                           -0.246 -0.354 -0.710 -0.748  0.416  0.575  0.702
##                             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
## m                            0.172
## s                            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