chapter 1 ros and max temp models
Sarah Kincy
2025-09-26
library and data:
library(ggplot2)
library(readr)
library(lme4)## Loading required package: Matrixlibrary(performance)
library(DHARMa)## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')library(see)
library(glmmTMB)
library(tweedie)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(sandwich) ##
## Attaching package: 'sandwich'## The following objects are masked from 'package:glmmTMB':
##
## meatHC, sandwichlibrary(lmtest) ## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(merDeriv) ## Loading required package: nonnest2## This is nonnest2 0.5-8.
## nonnest2 has not been tested with all combinations of supported model classes.## Loading required package: lavaan## This is lavaan 0.6-20
## lavaan is FREE software! Please report any bugs.library(clubSandwich)## Registered S3 methods overwritten by 'clubSandwich':
## method from
## bread.lmerMod merDeriv
## bread.mlm sandwichlibrary(emmeans)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'library(GLMMadaptive)##
## Attaching package: 'GLMMadaptive'## The following objects are masked from 'package:lavaan':
##
## anova, coef, fitted, nobs, predict, residuals## The following object is masked from 'package:lme4':
##
## negative.binomialchpt1_data <- read_csv("maindatachpt1.csv")## Rows: 51 Columns: 20## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): position, season, patch.location, herbicide, burn.date
## dbl (14): block, pair, max.temp, ros, live.biomass, litter.biomass, dead.bi...
## time (1): time stamp
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.chpt1_data <- chpt1_data %>%
mutate(
dead.biomass = ifelse(dead.biomass == 0, 0.0001, dead.biomass),
live.biomass = ifelse(live.biomass == 0, 0.0001, live.biomass)
)Global ROS model:
ros.global.model <- glmer(
ros ~ patch.location * season +
herbicide * season +
rel.humidity + mf.wind.sp +
live.biomass + litter.biomass + dead.biomass +
live.fmc + litter.fmc + dead.fmc +
(1 | block) + (1 | pair),
data = chpt1_data,
family = Gamma(link = "log"),
control = glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e5))
)
summary(ros.global.model)## Warning in vcov(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: ros ~ patch.location * season + herbicide * season + rel.humidity +
## mf.wind.sp + live.biomass + litter.biomass + dead.biomass +
## live.fmc + litter.fmc + dead.fmc + (1 | block) + (1 | pair)
## Data: chpt1_data
## Control: glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik -2*log(L) df.resid
## -719.3 -698.6 378.7 -757.3 3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.606e-04 -3.474e-05 4.530e-06 6.073e-05 1.184e-04
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 4.833e-01 6.952e-01
## block (Intercept) 2.790e-01 5.282e-01
## Residual 1.464e-09 3.826e-05
## Number of obs: 22, groups: pair, 11; block, 4
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) -2.3696757 1.3263122 -1.787 0.073992 .
## patch.locationinside 3.3579310 0.4844420 6.932 4.16e-12 ***
## patch.locationoutside 0.5109975 0.0443495 11.522 < 2e-16 ***
## seasongrowing 0.0000000 1.1665384 0.000 1.000000
## herbicideyes 0.3925705 0.5366362 0.732 0.464450
## rel.humidity -0.0892957 0.0231695 -3.854 0.000116 ***
## mf.wind.sp 1.4368704 0.4056300 3.542 0.000397 ***
## live.biomass 0.2717426 0.0297409 9.137 < 2e-16 ***
## litter.biomass -0.0455484 0.0096760 -4.707 2.51e-06 ***
## dead.biomass -0.1243885 0.0315565 -3.942 8.09e-05 ***
## live.fmc -0.0005734 0.0023690 -0.242 0.808742
## litter.fmc -0.0031709 0.0458279 -0.069 0.944837
## dead.fmc -0.0218253 0.0026284 -8.304 < 2e-16 ***
## patch.locationinside:seasongrowing -2.7807665 0.5938605 -4.683 2.83e-06 ***
## patch.locationoutside:seasongrowing -1.7553389 0.5054267 -3.473 0.000515 ***
## seasongrowing:herbicideyes -2.7786779 0.4772634 -5.822 5.81e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## optimizer (nloptwrap) convergence code: 0 (OK)
## Gradient contains NAsros_reduced <- update(ros.global.model, . ~ . - patch.location:season)
anova(ros_reduced, ros.global.model, test = "Chisq")## Data: chpt1_data
## Models:
## ros_reduced: ros ~ patch.location + season + herbicide + rel.humidity + mf.wind.sp + live.biomass + litter.biomass + dead.biomass + live.fmc + litter.fmc + dead.fmc + (1 | block) + (1 | pair) + season:herbicide
## ros.global.model: ros ~ patch.location * season + herbicide * season + rel.humidity + mf.wind.sp + live.biomass + litter.biomass + dead.biomass + live.fmc + litter.fmc + dead.fmc + (1 | block) + (1 | pair)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## ros_reduced 17 -692.52 -673.97 363.26 -726.52
## ros.global.model 19 -719.34 -698.61 378.67 -757.34 30.819 2 2.032e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ros_reduced2 <- update(ros.global.model, . ~ . - herbicide:season)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 2 negative eigenvaluesanova(ros_reduced2, ros.global.model, test = "Chisq")## Data: chpt1_data
## Models:
## ros_reduced2: ros ~ patch.location + season + herbicide + rel.humidity + mf.wind.sp + live.biomass + litter.biomass + dead.biomass + live.fmc + litter.fmc + dead.fmc + (1 | block) + (1 | pair) + patch.location:season
## ros.global.model: ros ~ patch.location * season + herbicide * season + rel.humidity + mf.wind.sp + live.biomass + litter.biomass + dead.biomass + live.fmc + litter.fmc + dead.fmc + (1 | block) + (1 | pair)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## ros_reduced2 18 -671.26 -651.62 353.63 -707.26
## ros.global.model 19 -719.34 -698.61 378.67 -757.34 50.076 1 1.479e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_res <- simulateResiduals(fittedModel = ros.global.model, n = 1000)## Simulations from your fitted model produce infinite values. Consider if this is sensible.plotQQunif(sim_res)
check_model(ros.global.model)## Simulations from your fitted model produce infinite values. Consider if this is sensible.
No fuel parameters:
ros.model1 <- glmer(
ros ~ patch.location * season +
herbicide * season +
rel.humidity + mf.wind.sp +
(1 | block) + (1 | pair),
data = chpt1_data,
family = Gamma(link = "log"),
control = glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e5))
)## boundary (singular) fit: see help('isSingular')summary(ros.model1)## Warning in vcov(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: ros ~ patch.location * season + herbicide * season + rel.humidity +
## mf.wind.sp + (1 | block) + (1 | pair)
## Data: chpt1_data
## Control: glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik -2*log(L) df.resid
## -241.9 -216.8 133.9 -267.9 38
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.22338 -0.61335 -0.07912 0.45094 2.72784
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 0.8032 0.8962
## block (Intercept) 0.0000 0.0000
## Residual 0.6049 0.7778
## Number of obs: 51, groups: pair, 20; block, 5
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) -5.79691 1.20498 -4.811 1.5e-06 ***
## patch.locationinside -0.31034 0.37337 -0.831 0.4059
## patch.locationoutside -0.24826 0.44102 -0.563 0.5735
## seasongrowing -1.33301 0.93726 -1.422 0.1550
## herbicideyes 0.35950 0.52182 0.689 0.4909
## rel.humidity 0.04109 0.03398 1.209 0.2266
## mf.wind.sp 0.23193 0.41362 0.561 0.5750
## patch.locationinside:seasongrowing 0.81264 0.67968 1.196 0.2318
## patch.locationoutside:seasongrowing -1.13080 0.73523 -1.538 0.1240
## seasongrowing:herbicideyes -2.22282 0.88189 -2.521 0.0117 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ptch.lctnn ptch.lctnt ssngrw hrbcdy rl.hmd mf.wn.
## ptch.lctnns -0.162
## ptch.lctnts -0.266 0.574
## seasongrwng 0.228 0.316 0.277
## herbicideys -0.226 0.130 0.537 0.161
## rel.humidty -0.597 -0.061 0.090 -0.595 0.133
## mf.wind.sp -0.137 0.022 -0.104 0.339 -0.128 -0.663
## ptch.lctnn: -0.003 -0.565 -0.281 -0.673 -0.025 0.333 -0.302
## ptch.lctnt: 0.094 -0.355 -0.577 -0.659 -0.291 0.150 -0.132
## ssngrwng:hr 0.053 -0.090 -0.289 -0.479 -0.552 0.177 -0.169
## ptch.lctnn: ptch.lctnt:
## ptch.lctnns
## ptch.lctnts
## seasongrwng
## herbicideys
## rel.humidty
## mf.wind.sp
## ptch.lctnn:
## ptch.lctnt: 0.665
## ssngrwng:hr 0.239 0.520
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')ros_reduced <- update(ros.model1, . ~ . - patch.location:season)## boundary (singular) fit: see help('isSingular')anova(ros_reduced, ros.model1, test = "Chisq")## Data: chpt1_data
## Models:
## ros_reduced: ros ~ patch.location + season + herbicide + rel.humidity + mf.wind.sp + (1 | block) + (1 | pair) + season:herbicide
## ros.model1: ros ~ patch.location * season + herbicide * season + rel.humidity + mf.wind.sp + (1 | block) + (1 | pair)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## ros_reduced 11 -237.75 -216.50 129.88 -259.75
## ros.model1 13 -241.87 -216.76 133.94 -267.87 8.1179 2 0.01727 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ros_reduced2 <- update(ros.model1, . ~ . - herbicide:season)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluesanova(ros_reduced2, ros.model1, test = "Chisq")## Data: chpt1_data
## Models:
## ros_reduced2: ros ~ patch.location + season + herbicide + rel.humidity + mf.wind.sp + (1 | block) + (1 | pair) + patch.location:season
## ros.model1: ros ~ patch.location * season + herbicide * season + rel.humidity + mf.wind.sp + (1 | block) + (1 | pair)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## ros_reduced2 12 -238.17 -214.98 131.08 -262.17
## ros.model1 13 -241.87 -216.76 133.94 -267.87 5.705 1 0.01692 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_res <- simulateResiduals(fittedModel = ros.model1, n = 1000)
plotQQunif(sim_res)
check_model(ros.model1)
Removing season * location and putting fuel parameters back in.
chpt1_data_scaled <- chpt1_data
numeric_vars <- c("rel.humidity", "mf.wind.sp",
"live.biomass", "litter.biomass", "dead.biomass",
"live.fmc", "litter.fmc", "dead.fmc")
chpt1_data_scaled[numeric_vars] <- scale(chpt1_data_scaled[numeric_vars])
ros.model2 <- glmer(
ros ~ patch.location +
herbicide * season +
rel.humidity + mf.wind.sp +
live.biomass + litter.biomass + dead.biomass +
live.fmc + litter.fmc + dead.fmc +
(1 | block) + (1 | pair),
data = chpt1_data_scaled,
family = Gamma(link = "log"),
control = glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e5))
)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Problem with Hessian check (infinite or missing values?)summary(ros.model2)## Warning in vcov(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula:
## ros ~ patch.location + herbicide * season + rel.humidity + mf.wind.sp +
## live.biomass + litter.biomass + dead.biomass + live.fmc +
## litter.fmc + dead.fmc + (1 | block) + (1 | pair)
## Data: chpt1_data_scaled
## Control: glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik -2*log(L) df.resid
## -680.7 -662.1 357.3 -714.7 5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.717e-04 -5.184e-05 7.340e-06 6.243e-05 3.793e-04
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 2.152e+00 1.467e+00
## block (Intercept) 3.469e-03 5.890e-02
## Residual 1.178e-09 3.432e-05
## Number of obs: 22, groups: pair, 11; block, 4
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) -3.651955 0.614507 -5.943 2.80e-09 ***
## patch.locationinside 3.628739 0.007079 512.616 < 2e-16 ***
## patch.locationoutside 0.234258 0.025143 9.317 < 2e-16 ***
## herbicideyes -3.703189 0.116567 -31.769 < 2e-16 ***
## seasongrowing -4.171922 1.066120 -3.913 9.11e-05 ***
## rel.humidity 0.429089 0.145933 2.940 0.00328 **
## mf.wind.sp -0.673620 0.285347 -2.361 0.01824 *
## live.biomass 0.921095 0.003144 292.956 < 2e-16 ***
## litter.biomass -1.943671 0.138451 -14.039 < 2e-16 ***
## dead.biomass 0.199462 0.098082 2.034 0.04199 *
## live.fmc 0.882916 0.036466 24.212 < 2e-16 ***
## litter.fmc 3.022834 0.031468 96.062 < 2e-16 ***
## dead.fmc -0.467210 0.038796 -12.043 < 2e-16 ***
## herbicideyes:seasongrowing -0.757294 0.119649 -6.329 2.46e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Problem with Hessian check (infinite or missing values?)ros_reduced2 <- update(ros.model2, . ~ . - herbicide:season)## boundary (singular) fit: see help('isSingular')anova(ros_reduced2, ros.model2, test = "Chisq")## Data: chpt1_data_scaled
## Models:
## ros_reduced2: ros ~ patch.location + herbicide + season + rel.humidity + mf.wind.sp + live.biomass + litter.biomass + dead.biomass + live.fmc + litter.fmc + dead.fmc + (1 | block) + (1 | pair)
## ros.model2: ros ~ patch.location + herbicide * season + rel.humidity + mf.wind.sp + live.biomass + litter.biomass + dead.biomass + live.fmc + litter.fmc + dead.fmc + (1 | block) + (1 | pair)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## ros_reduced2 16 -129.97 -112.51 80.99 -161.97
## ros.model2 17 -680.68 -662.13 357.34 -714.68 552.71 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_res <- simulateResiduals(fittedModel = ros.model2, n = 1000)## Simulations from your fitted model produce infinite values. Consider if this is sensible.plotQQunif(sim_res)
check_model(ros.model2)## Simulations from your fitted model produce infinite values. Consider if this is sensible.
Removing block random effect:
chpt1_data_scaled <- chpt1_data
numeric_vars <- c("rel.humidity", "mf.wind.sp",
"live.biomass", "litter.biomass", "dead.biomass",
"live.fmc", "litter.fmc", "dead.fmc")
chpt1_data_scaled[numeric_vars] <- scale(chpt1_data_scaled[numeric_vars])
ros.model3 <- glmer(
ros ~ patch.location +
herbicide * season +
rel.humidity + mf.wind.sp +
live.biomass + litter.biomass + dead.biomass +
live.fmc + litter.fmc + dead.fmc +
(1 | pair),
data = chpt1_data_scaled,
family = Gamma(link = "log"),
control = glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e5))
)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluessummary(ros.model3)## Warning in vcov(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula:
## ros ~ patch.location + herbicide * season + rel.humidity + mf.wind.sp +
## live.biomass + litter.biomass + dead.biomass + live.fmc +
## litter.fmc + dead.fmc + (1 | pair)
## Data: chpt1_data_scaled
## Control: glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik -2*log(L) df.resid
## -673.8 -656.3 352.9 -705.8 6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.714e-04 -4.592e-05 0.000e+00 4.128e-05 2.770e-04
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 5.309e+00 2.304e+00
## Residual 3.635e-09 6.029e-05
## Number of obs: 22, groups: pair, 11
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) -3.655632 0.962884 -3.797 0.000147 ***
## patch.locationinside 3.628826 0.011095 327.064 < 2e-16 ***
## patch.locationoutside 0.233948 0.039407 5.937 2.91e-09 ***
## herbicideyes -3.704621 0.182696 -20.278 < 2e-16 ***
## seasongrowing 0.000000 1.669681 0.000 1.000000
## rel.humidity 0.427296 0.228720 1.868 0.061734 .
## mf.wind.sp -0.677127 0.447224 -1.514 0.130009
## live.biomass 0.921134 0.004928 186.923 < 2e-16 ***
## litter.biomass -1.941969 0.216995 -8.949 < 2e-16 ***
## dead.biomass 0.200668 0.153725 1.305 0.191767
## live.fmc 0.883364 0.057153 15.456 < 2e-16 ***
## litter.fmc 3.022448 0.049319 61.283 < 2e-16 ***
## dead.fmc -0.467686 0.060805 -7.692 1.45e-14 ***
## herbicideyes:seasongrowing -0.755823 0.187526 -4.030 5.57e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluesros_reduced2 <- update(ros.model3, . ~ . - herbicide:season)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 1.86015 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?anova(ros_reduced2, ros.model3, test = "Chisq")## Data: chpt1_data_scaled
## Models:
## ros_reduced2: ros ~ patch.location + herbicide + season + rel.humidity + mf.wind.sp + live.biomass + litter.biomass + dead.biomass + live.fmc + litter.fmc + dead.fmc + (1 | pair)
## ros.model3: ros ~ patch.location + herbicide * season + rel.humidity + mf.wind.sp + live.biomass + litter.biomass + dead.biomass + live.fmc + litter.fmc + dead.fmc + (1 | pair)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## ros_reduced2 15 -292.27 -275.90 161.13 -322.27
## ros.model3 16 -673.77 -656.31 352.88 -705.77 383.5 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_res <- simulateResiduals(fittedModel = ros.model3, n = 1000)## Simulations from your fitted model produce infinite values. Consider if this is sensible.plotQQunif(sim_res)
check_model(ros.model3)## Simulations from your fitted model produce infinite values. Consider if this is sensible.
chpt1_data_scaled <- chpt1_data
numeric_vars <- c("rel.humidity", "mf.wind.sp",
"live.biomass", "litter.biomass", "dead.biomass",
"live.fmc", "litter.fmc", "dead.fmc")
chpt1_data_scaled[numeric_vars] <- scale(chpt1_data_scaled[numeric_vars])
ros.model4 <- glmer(
ros ~ patch.location +
herbicide * season +
rel.humidity + mf.wind.sp +
dead.biomass + dead.fmc +
(1 | pair),
data = chpt1_data_scaled,
family = Gamma(link = "log"),
control = glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e5))
)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00494702 (tol = 0.002, component 1)summary(ros.model4)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula:
## ros ~ patch.location + herbicide * season + rel.humidity + mf.wind.sp +
## dead.biomass + dead.fmc + (1 | pair)
## Data: chpt1_data_scaled
## Control: glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik -2*log(L) df.resid
## -153.5 -136.7 88.8 -177.5 18
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.28511 -0.41812 -0.04706 0.48194 1.82801
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 0.5384 0.7337
## Residual 0.2286 0.4781
## Number of obs: 30, groups: pair, 14
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) -3.6936 0.4390 -8.413 < 2e-16 ***
## patch.locationinside 0.9544 0.2308 4.135 3.55e-05 ***
## patch.locationoutside -0.2699 0.2850 -0.947 0.343666
## herbicideyes 0.4177 0.4368 0.956 0.338885
## seasongrowing -1.0761 0.8300 -1.297 0.194775
## rel.humidity 0.4064 0.3043 1.336 0.181698
## mf.wind.sp 0.1662 0.2472 0.672 0.501333
## dead.biomass -0.8344 0.1658 -5.033 4.84e-07 ***
## dead.fmc -0.5312 0.1170 -4.542 5.58e-06 ***
## herbicideyes:seasongrowing -2.0633 0.5499 -3.752 0.000175 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ptch.lctnn ptch.lctnt hrbcdy ssngrw rl.hmd mf.wn. dd.bms
## ptch.lctnns -0.343
## ptch.lctnts -0.380 0.320
## herbicideys -0.347 0.135 0.597
## seasongrwng -0.511 0.029 -0.066 -0.044
## rel.humidty 0.236 0.017 0.121 0.111 -0.511
## mf.wind.sp 0.004 -0.254 -0.168 -0.067 0.123 -0.216
## dead.biomss 0.228 -0.508 -0.160 -0.499 0.033 -0.009 0.047
## dead.fmc 0.188 -0.287 0.064 0.117 -0.274 0.277 -0.079 0.279
## hrbcdys:ssn 0.285 -0.251 0.027 -0.412 -0.303 0.320 -0.261 0.448
## dd.fmc
## ptch.lctnns
## ptch.lctnts
## herbicideys
## seasongrwng
## rel.humidty
## mf.wind.sp
## dead.biomss
## dead.fmc
## hrbcdys:ssn 0.392
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00494702 (tol = 0.002, component 1)ros_reduced2 <- update(ros.model4, . ~ . - herbicide:season)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0297959 (tol = 0.002, component 1)anova(ros_reduced2, ros.model4, test = "Chisq")## Data: chpt1_data_scaled
## Models:
## ros_reduced2: ros ~ patch.location + herbicide + season + rel.humidity + mf.wind.sp + dead.biomass + dead.fmc + (1 | pair)
## ros.model4: ros ~ patch.location + herbicide * season + rel.humidity + mf.wind.sp + dead.biomass + dead.fmc + (1 | pair)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## ros_reduced2 11 -147.20 -131.79 84.602 -169.20
## ros.model4 12 -153.51 -136.70 88.757 -177.51 8.3105 1 0.003942 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_res <- simulateResiduals(fittedModel = ros.model4, n = 1000)
plotQQunif(sim_res)
check_model(ros.model4)