Module 5 Content Review
Camie Tomlinson
2022-10-09
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(lme4))##Load in Data
schooldata <- haven::read_dta("schoolsmoke-1.dta")
glimpse(schooldata)## Rows: 1,600
## Columns: 5
## $ school <dbl> 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 19…
## $ thk <dbl> 1, 2, 3, 1, 2, 2, 3, 3, 2, 2, 4, 2, 2, 3, 4, 4, 1, 3, 4, 2, 1,…
## $ prethk <dbl> 5, 3, 4, 0, 1, 3, 1, 2, 3, 2, 2, 2, 3, 3, 2, 3, 3, 3, 2, 1, 3,…
## $ cc <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ gprethk <dbl> 2.153846, 2.153846, 2.153846, 2.153846, 2.153846, 2.153846, 2.…#Changing cc to a factor variable/telling R it is categorical
schooldata.1 <- schooldata %>%
mutate(., cc.factor = as_factor(cc))
glimpse(schooldata.1) ## Rows: 1,600
## Columns: 6
## $ school <dbl> 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, 193, …
## $ thk <dbl> 1, 2, 3, 1, 2, 2, 3, 3, 2, 2, 4, 2, 2, 3, 4, 4, 1, 3, 4, 2, …
## $ prethk <dbl> 5, 3, 4, 0, 1, 3, 1, 2, 3, 2, 2, 2, 3, 3, 2, 3, 3, 3, 2, 1, …
## $ cc <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ gprethk <dbl> 2.153846, 2.153846, 2.153846, 2.153846, 2.153846, 2.153846, …
## $ cc.factor <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …##Null Model
model.null <- lmer(thk ~ (1|school), REML = FALSE, data = schooldata.1)
summary(model.null)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: thk ~ (1 | school)
## Data: schooldata.1
##
## AIC BIC logLik deviance df.resid
## 4833.0 4849.1 -2413.5 4827.0 1597
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9039 -0.7734 0.0693 0.9240 1.7430
##
## Random effects:
## Groups Name Variance Std.Dev.
## school (Intercept) 0.09438 0.3072
## Residual 1.16219 1.0780
## Number of obs: 1600, groups: school, 28
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.60750 0.06538 39.88#Calculate ICC
null.ICC <- 0.09438 /(0.09438 + 1.16219 )
null.ICC## [1] 0.07510923#turning off scientific notation
options(scipen=999)lmerTest::rand(model.null)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## thk ~ (1 | school)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 3 -2413.5 4833.0
## (1 | school) 2 -2445.6 4895.2 64.16 1 0.000000000000001147 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##Conditional Random Intercept Models
model.1 <- lmer(thk ~ prethk + (1|school), REML = FALSE, data = schooldata.1)
summary(model.1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: thk ~ prethk + (1 | school)
## Data: schooldata.1
##
## AIC BIC logLik deviance df.resid
## 4733.5 4755.0 -2362.7 4725.5 1596
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.12273 -0.82727 0.03082 0.83756 2.15165
##
## Random effects:
## Groups Name Variance Std.Dev.
## school (Intercept) 0.07519 0.2742
## Residual 1.09321 1.0456
## Number of obs: 1600, groups: school, 28
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.15164 0.07414 29.02
## prethk 0.21747 0.02122 10.25
##
## Correlation of Fixed Effects:
## (Intr)
## prethk -0.598model.2 <- lmer(thk ~ prethk + cc.factor + (1|school), REML = FALSE, data = schooldata.1)
summary(model.2)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: thk ~ prethk + cc.factor + (1 | school)
## Data: schooldata.1
##
## AIC BIC logLik deviance df.resid
## 4721.3 4748.2 -2355.6 4711.3 1595
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.20705 -0.80300 0.04029 0.80179 2.15135
##
## Random effects:
## Groups Name Variance Std.Dev.
## school (Intercept) 0.03333 0.1826
## Residual 1.09412 1.0460
## Number of obs: 1600, groups: school, 28
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.94507 0.07699 25.265
## prethk 0.22262 0.02113 10.535
## cc.factor1 0.39215 0.08944 4.385
##
## Correlation of Fixed Effects:
## (Intr) prethk
## prethk -0.585
## cc.factor1 -0.580 0.023model.3 <- lmer(thk ~ prethk + cc.factor + gprethk + (1|school), REML = FALSE, data = schooldata.1)
summary(model.3)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: thk ~ prethk + cc.factor + gprethk + (1 | school)
## Data: schooldata.1
##
## AIC BIC logLik deviance df.resid
## 4712.6 4744.9 -2350.3 4700.6 1594
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.23332 -0.79567 0.02374 0.79254 2.16704
##
## Random effects:
## Groups Name Variance Std.Dev.
## school (Intercept) 0.01734 0.1317
## Residual 1.09334 1.0456
## Number of obs: 1600, groups: school, 28
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.07609 0.25478 4.224
## prethk 0.21260 0.02138 9.944
## cc.factor1 0.43075 0.07554 5.702
## gprethk 0.41835 0.11930 3.507
##
## Correlation of Fixed Effects:
## (Intr) prethk cc.fc1
## prethk 0.000
## cc.factor1 -0.288 0.000
## gprethk -0.963 -0.179 0.148#Summarizing conditional random intercept model
library(modelsummary)
library(broom.mixed)
models <- list(model.null, model.1, model.2, model.3)
modelsummary(models, output = "markdown")| Model 1 | Model 2 | Model 3 | Model 4 | |
|---|---|---|---|---|
| (Intercept) | 2.607 | 2.152 | 1.945 | 1.076 |
| (0.065) | (0.074) | (0.077) | (0.255) | |
| prethk | 0.217 | 0.223 | 0.213 | |
| (0.021) | (0.021) | (0.021) | ||
| cc.factor1 | 0.392 | 0.431 | ||
| (0.089) | (0.076) | |||
| gprethk | 0.418 | |||
| (0.119) | ||||
| SD (Intercept school) | 0.307 | 0.274 | 0.183 | 0.132 |
| SD (Observations) | 1.078 | 1.046 | 1.046 | 1.046 |
| Num.Obs. | 1600 | 1600 | 1600 | 1600 |
| R2 Marg. | 0.000 | 0.060 | 0.091 | 0.109 |
| R2 Cond. | 0.075 | 0.121 | 0.118 | 0.123 |
| AIC | 4836.6 | 4743.2 | 4734.5 | 4729.0 |
| BIC | 4852.7 | 4764.7 | 4761.4 | 4761.3 |
| ICC | 0.1 | 0.1 | 0.0 | 0.0 |
| RMSE | 1.07 | 1.04 | 1.04 | 1.04 |
##Adding Random Slope - PRETHK at level 2
model.3.prethk_slp <- lmer(thk ~ prethk + cc.factor + gprethk + (prethk|school), REML = FALSE, data = schooldata.1)
summary(model.3.prethk_slp)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: thk ~ prethk + cc.factor + gprethk + (prethk | school)
## Data: schooldata.1
##
## AIC BIC logLik deviance df.resid
## 4715.5 4758.6 -2349.8 4699.5 1592
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.32829 -0.81670 0.02938 0.81294 2.14320
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## school (Intercept) 0.049006 0.22137
## prethk 0.002063 0.04542 -0.97
## Residual 1.090740 1.04438
## Number of obs: 1600, groups: school, 28
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.90319 0.24823 3.639
## prethk 0.21136 0.02323 9.097
## cc.factor1 0.42343 0.07337 5.771
## gprethk 0.50405 0.11449 4.403
##
## Correlation of Fixed Effects:
## (Intr) prethk cc.fc1
## prethk -0.065
## cc.factor1 -0.295 -0.004
## gprethk -0.951 -0.174 0.160lmerTest::rand(model.3.prethk_slp)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## thk ~ prethk + cc.factor + gprethk + (prethk | school)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 -2349.8 4715.5
## prethk in (prethk | school) 6 -2350.3 4712.6 1.0762 2 0.5839##Testing a cross-level interaction
#centering predictors
schooldata.2 <- schooldata.1 %>%
mutate(.,
prethk_MC = prethk - mean(prethk),
gprethk_MC = gprethk - mean(gprethk))
psych::describe(schooldata.2, fast = TRUE)## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf## vars n mean sd min max range se
## school 1 1600 421.94 112.67 193.00 515.00 322.00 2.82
## thk 2 1600 2.59 1.12 1.00 4.00 3.00 0.03
## prethk 3 1600 2.07 1.26 0.00 6.00 6.00 0.03
## cc 4 1600 0.48 0.50 0.00 1.00 1.00 0.01
## gprethk 5 1600 2.07 0.30 1.59 3.06 1.47 0.01
## cc.factor 6 1600 NaN NA Inf -Inf -Inf NA
## prethk_MC 7 1600 0.00 1.26 -2.07 3.93 6.00 0.03
## gprethk_MC 8 1600 0.00 0.30 -0.48 0.99 1.47 0.01#run main effects model
model.4 <- lmer(thk ~ prethk_MC + cc.factor + gprethk_MC + (1|school), REML = FALSE, data = schooldata.2)
summary(model.4)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: thk ~ prethk_MC + cc.factor + gprethk_MC + (1 | school)
## Data: schooldata.2
##
## AIC BIC logLik deviance df.resid
## 4712.6 4744.9 -2350.3 4700.6 1594
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.23332 -0.79567 0.02374 0.79254 2.16704
##
## Random effects:
## Groups Name Variance Std.Dev.
## school (Intercept) 0.01734 0.1317
## Residual 1.09334 1.0456
## Number of obs: 1600, groups: school, 28
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.38176 0.05235 45.499
## prethk_MC 0.21260 0.02138 9.944
## cc.factor1 0.43075 0.07554 5.702
## gprethk_MC 0.41835 0.11930 3.507
##
## Correlation of Fixed Effects:
## (Intr) prt_MC cc.fc1
## prethk_MC 0.000
## cc.factor1 -0.701 0.000
## gprethk_MC -0.123 -0.179 0.148#interaction effects model
model.5.intax <- lmer(thk ~ prethk_MC + cc.factor + gprethk_MC + prethk_MC:gprethk_MC + (1|school), REML = FALSE, data = schooldata.2)
summary(model.5.intax)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: thk ~ prethk_MC + cc.factor + gprethk_MC + prethk_MC:gprethk_MC +
## (1 | school)
## Data: schooldata.2
##
## AIC BIC logLik deviance df.resid
## 4702.9 4740.5 -2344.4 4688.9 1593
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.24644 -0.81395 0.04367 0.80525 2.07026
##
## Random effects:
## Groups Name Variance Std.Dev.
## school (Intercept) 0.0165 0.1284
## Residual 1.0857 1.0420
## Number of obs: 1600, groups: school, 28
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.35901 0.05205 45.322
## prethk_MC 0.20797 0.02135 9.743
## cc.factor1 0.43001 0.07453 5.770
## gprethk_MC 0.36513 0.11890 3.071
## prethk_MC:gprethk_MC 0.22898 0.06663 3.437
##
## Correlation of Fixed Effects:
## (Intr) prt_MC cc.fc1 gpr_MC
## prethk_MC 0.008
## cc.factor1 -0.695 0.000
## gprethk_MC -0.104 -0.171 0.148
## prth_MC:_MC -0.128 -0.063 -0.003 -0.130#chi-square diff test – include interaction effect?
anova(model.4, model.5.intax)## Data: schooldata.2
## Models:
## model.4: thk ~ prethk_MC + cc.factor + gprethk_MC + (1 | school)
## model.5.intax: thk ~ prethk_MC + cc.factor + gprethk_MC + prethk_MC:gprethk_MC + (1 | school)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.4 6 4712.6 4744.9 -2350.3 4700.6
## model.5.intax 7 4702.9 4740.5 -2344.4 4688.9 11.764 1 0.000604 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#visualize moderation
interplot::interplot(model.5.intax, var1 = "gprethk_MC", var2 = "prethk_MC")
library(modelsummary)
library(broom.mixed)
models.2 <- list(model.4, model.5.intax)
modelsummary(models.2, output = "markdown")| Model 1 | Model 2 | |
|---|---|---|
| (Intercept) | 2.382 | 2.359 |
| (0.052) | (0.052) | |
| prethk_MC | 0.213 | 0.208 |
| (0.021) | (0.021) | |
| cc.factor1 | 0.431 | 0.430 |
| (0.076) | (0.075) | |
| gprethk_MC | 0.418 | 0.365 |
| (0.119) | (0.119) | |
| prethk_MCgprethk_MC | 0.229 | |
| (0.067) | ||
| SD (Intercept school) | 0.132 | 0.128 |
| SD (Observations) | 1.046 | 1.042 |
| Num.Obs. | 1600 | 1600 |
| R2 Marg. | 0.109 | 0.115 |
| R2 Cond. | 0.123 | 0.129 |
| AIC | 4729.0 | 4722.9 |
| BIC | 4761.3 | 4760.5 |
| ICC | 0.0 | 0.0 |
| RMSE | 1.04 | 1.04 |
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.