Baseline model
Intercept-only model to estimate ICC — how much of the variation in
darkness is due to dojo differences.
baseline = lmer(darkness ~ 1 + (1 | dojo_id), data = d)
summary(baseline)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: darkness ~ 1 + (1 | dojo_id)
#> Data: d
#>
#> REML criterion at convergence: 10349.2
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -3.0950 -0.6789 -0.0096 0.6590 3.7008
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> dojo_id (Intercept) 0.700 0.8367
#> Residual 3.261 1.8059
#> Number of obs: 2530, groups: dojo_id, 95
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 5.00835 0.09389 53.34
Fixed anger effect
Add anger as a fixed effect to test whether anger predicts darkness
across all dojos.
fixed_anger = lmer(darkness ~ 1 + anger + (1 | dojo_id), data = d)
summary(fixed_anger)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: darkness ~ 1 + anger + (1 | dojo_id)
#> Data: d
#>
#> REML criterion at convergence: 9150.9
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -3.7574 -0.6651 -0.0064 0.6657 4.5388
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> dojo_id (Intercept) 0.7545 0.8686
#> Residual 1.9883 1.4101
#> Number of obs: 2530, groups: dojo_id, 95
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 0.484524 0.148227 3.269
#> anger 0.377712 0.009564 39.495
#>
#> Correlation of Fixed Effects:
#> (Intr)
#> anger -0.773
Random + fixed anger effect
Allow the effect of anger to vary across dojos (random slope).
fixed_random_anger = lmer(darkness ~ 1 + anger + (1 + anger | dojo_id), data = d)
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00757766 (tol = 0.002, component 1)
#> See ?lme4::convergence and ?lme4::troubleshooting.
summary(fixed_random_anger)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: darkness ~ 1 + anger + (1 + anger | dojo_id)
#> Data: d
#>
#> REML criterion at convergence: 9092.9
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -3.9545 -0.6485 -0.0211 0.6749 4.7032
#>
#> Random effects:
#> Groups Name Variance Std.Dev. Corr
#> dojo_id (Intercept) 2.63690 1.6239
#> anger 0.01252 0.1119 -0.84
#> Residual 1.87626 1.3698
#> Number of obs: 2530, groups: dojo_id, 95
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 0.45873 0.20698 2.216
#> anger 0.37912 0.01518 24.977
#>
#> Correlation of Fixed Effects:
#> (Intr)
#> anger -0.891
#> optimizer (nloptwrap) convergence code: 0 (OK)
#> Model failed to converge with max|grad| = 0.00757766 (tol = 0.002, component 1)
#> See ?lme4::convergence and ?lme4::troubleshooting.
Main model: does emotional attachment increase darkness, controlling
for anger?
model = lmer(darkness ~ anger + emotional_bonds +
(anger + emotional_bonds | dojo_id), data = d,
control = lmerControl(optimizer = "bobyqa"))
summary(model)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: darkness ~ anger + emotional_bonds + (anger + emotional_bonds |
#> dojo_id)
#> Data: d
#> Control: lmerControl(optimizer = "bobyqa")
#>
#> REML criterion at convergence: 8789.8
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -4.4330 -0.6266 0.0003 0.6543 3.8001
#>
#> Random effects:
#> Groups Name Variance Std.Dev. Corr
#> dojo_id (Intercept) 6.91008 2.6287
#> anger 0.01302 0.1141 -0.62
#> emotional_bonds 0.02215 0.1488 -0.78 0.11
#> Residual 1.56219 1.2499
#> Number of obs: 2530, groups: dojo_id, 95
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) -0.84430 0.31363 -2.692
#> anger 0.37820 0.01492 25.343
#> emotional_bonds 0.10897 0.01791 6.083
#>
#> Correlation of Fixed Effects:
#> (Intr) anger
#> anger -0.638
#> emotnl_bnds -0.756 0.071
estimates(model)
#> refitting model(s) with ML (instead of REML)
#> Fixed Effects:
#> (Intercept) anger emotional_bonds
#> -0.8443012 0.3782025 0.1089741
#>
#>
#> Random Effects:
#> Groups Name Std.Dev. Corr
#> dojo_id (Intercept) 2.62870
#> anger 0.11413 -0.623
#> emotional_bonds 0.14882 -0.777 0.111
#> Residual 1.24987
#>
#>
#> ICC and Design Effect:
#> icc design.effect
#> 0.1767232 5.5296945
#>
#>
#> R Squared:
#>
#> (Intercept) Residual
#> -8.8710696 0.5209714
Visualise model predictions (overall)
Shows the overall model-predicted relationship between predictors and
darkness.
visualize(model, plot = "model")
#> Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
#> ℹ Please use tidy evaluation idioms with `aes()`.
#> ℹ See also `vignette("ggplot2-in-packages")` for more information.
#> ℹ The deprecated feature was likely used in the flexplot package.
#> Please report the issue to the authors.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
Visualise emotional bonds effect by anger level
Shows how emotional_bonds predicts darkness at different levels of
anger, across dojos.
visualize(model, plot = "model",
formula = darkness ~ emotional_bonds + dojo_id | anger)
Visualise emotional bonds across 30 sampled dojos
Plots the emotional_bonds effect for a random sample of 30 dojos to
avoid overplotting.
visualize(model, plot = "model",
formula = darkness ~ emotional_bonds + dojo_id, sample = 30)
#> Note: You didn't choose to plot anger so I am inputting the mean
#> It looks like you're trying to plot more than 6 colors/lines/symbols.
#> I gotta give it to you...you're ambitious. Alas, I can't do that, so I'm removing the colors/lines/symbols.
#> I hope we can still be friends.
Visualise anger across 30 sampled dojos
Plots the anger effect for a random sample of 30 dojos.
visualize(model, plot = "model",
formula = darkness ~ anger + dojo_id, sample = 30)
#> Note: You didn't choose to plot emotional_bonds so I am inputting the mean
#> It looks like you're trying to plot more than 6 colors/lines/symbols.
#> I gotta give it to you...you're ambitious. Alas, I can't do that, so I'm removing the colors/lines/symbols.
#> I hope we can still be friends.
Model comparison: is the random slope for anger necessary?
Compare full model (random slopes for both) vs reduced model (random
slope for emotional_bonds only).
reduced = lmer(darkness ~ emotional_bonds + anger +
(emotional_bonds | dojo_id), data = d,
control = lmerControl(optimizer = "bobyqa"))
compare.fits(darkness ~ anger | emotional_bonds + dojo_id,
data = d, model1 = model, model2 = reduced, re = T)
model.comparison(model, reduced)
#> refitting model(s) with ML (instead of REML)
#> $statistics
#> aic bic bayes.factor p
#> model 8809.836 8868.196 118335836782 <2e-16
#> reduced 8878.337 8919.189 0
#>
#> $predicted_differences
#> 0% 25% 50% 75% 100%
#> 0.000 0.028 0.091 0.230 2.062
#>
#> $r_squared_change
#> (Intercept) emotional_bonds Residual
#> -0.58638018 0.01003077 0.06820976
