Module 10
Caitlin Mayton
This week we used ProjectSTAR data again, which is…“a sample of more than 6,000 students nested within 200+ schools collected between 1985 and 1990 as part of the Tennessee Class Size Experiment.” -Dr. Broda (Module 7)
#Load in our MVP packages
suppressPackageStartupMessages(library(tidyverse))
library(lme4)
library(Hmisc)Part 1: loading in the data
projectSTAR <- haven::read_dta("projectSTAR.dta")
glimpse(projectSTAR)Rows: 6,325
Columns: 27
$ stdntid <dbl> 10001, 10133, 10246, 10263, 10266, 10275, 10281, 10282, 10285, 10286, 10287, 10289, 10290, 10291, 10292, 10294, 10295, 10299, 10300, 10302, 1…
$ race <dbl+lbl> 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2,…
$ gender <dbl+lbl> 1, 1, 2, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 1,…
$ FLAGSGK <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ flaggk <dbl+lbl> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,…
$ gkclasstype <dbl+lbl> 3, 3, 3, 1, 2, 3, 1, 3, 1, 2, 3, 3, 2, 3, 3, 3, 1, 2, 3, 3, 2, 1, 1, 3, 2, 1, 1, 2, 3, 2, 3, 3, 3, 2, 3, 2, 3, 3, 1, 1, 3, 1, 3, 1, 3, 1,…
$ gkschid <dbl> 169229, 169280, 218562, 205492, 257899, 161176, 189382, 189382, 201449, 230612, 128068, 262937, 169219, 169229, 205491, 159171, 161176, 16821…
$ gksurban <dbl+lbl> 2, 2, 4, 2, 3, 3, 2, 2, 3, 3, 3, 2, 2, 2, 2, 3, 3, 3, 3, 3, 2, 1, 2, 3, 1, 4, 2, 3, 2, 3, 4, 3, 2, 2, 2, 3, 2, 2, 4, 3, 1, 1, 1, 1, 1, 1,…
$ gktchid <dbl> 16922904, 16928003, 21856202, 20549204, 25789904, 16117602, 18938204, 18938203, 20144901, 23061203, 12806803, 26293701, 16921902, 16922904, 2…
$ gktgen <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
$ gktrace <dbl+lbl> 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, NA, 1, 1, 1, …
$ gkthighdegree <dbl+lbl> 2, 2, 3, 2, 3, 2, 3, 3, 3, 3, 2, 2, 3, 2, 2, 2, 2, 2, 2, 3, 2, 2, 3, 2, 2, 3, 3, 2, 2, 2, NA, 2, 2, 3, …
$ gktcareer <dbl+lbl> 4, 4, 4, 4, 4, NA, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, NA, 2, 4, 4, 2, 4, 4, 4, 4, NA, NA, 4, 4, 6, …
$ gktyears <dbl> 5, 7, 8, 3, 12, 2, 7, 14, 4, 6, 11, 16, 12, 5, 17, 10, 6, 10, 13, 9, 18, 1, 7, 13, 1, 14, 18, 8, 8, 4, NA, 13, 17, 15, 16, 13, 3, 17, 11, 2, …
$ gkclasssize <dbl> 24, 22, 17, 17, 24, 24, 13, 24, 14, 23, 23, 22, 20, 24, 23, 27, 17, 24, 22, 23, 22, 17, 13, 18, 20, 14, 14, 23, 24, 24, 21, 18, 23, 19, 22, 2…
$ gkfreelunch <dbl+lbl> 2, 2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 2, 1, 2, 2, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1,…
$ gkrepeat <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ gkspeced <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 1, 2, 2,…
$ gkspecin <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
$ gkpresent <dbl> 161, 175, 132, 178, 170, 94, 160, 154, 172, 95, 163, 172, 180, 149, 173, NA, 160, 173, 170, 161, 128, 151, 64, 84, 145, 160, 161, 177, 154, 1…
$ gkabsent <dbl> 19, 5, 28, 2, 10, 3, 2, 7, 8, 2, 17, 8, 0, 31, 7, NA, 20, 7, 2, 19, 24, 0, 8, 9, 7, 4, 19, 3, 3, 4, 15, 26, 6, 5, 10, 3, 15, 13, 2, 1, 12, 7,…
$ gktreadss <dbl> NA, 427, 450, 483, 456, 411, 443, 448, 463, 472, 428, 545, 408, 422, 472, NA, 437, 450, 497, 460, 408, NA, 467, 403, 411, 490, 503, NA, 442, …
$ gktmathss <dbl> NA, 478, 494, 513, 513, 468, 473, 449, 520, 536, 484, 626, 454, 439, 528, NA, 463, 484, 547, 576, 449, NA, 576, 392, 444, 559, 576, NA, 547, …
$ gktlistss <dbl> NA, 509, 549, 554, 520, 571, 595, 540, 565, 595, NA, 622, 474, 536, 578, NA, 554, 565, 524, 595, 512, NA, 565, 501, 505, 595, 671, NA, 532, 5…
$ gkwordskillss <dbl> NA, 418, 444, 431, 468, 396, 444, 444, 480, 486, 423, 524, 410, 423, 458, NA, 427, 435, 512, 449, 405, NA, 480, 410, 405, 512, 480, NA, 435, …
$ gkmotivraw <dbl> 23, 24, 28, 27, 25, 24, NA, NA, 26, 27, 24, 24, 23, 28, 24, NA, 26, 25, NA, 26, 25, 29, NA, 24, 26, 24, 25, 28, NA, 27, NA, 17, 24, 23, 25, 2…
$ gkselfconcraw <dbl> 52, 53, 56, 61, 54, 55, NA, NA, 52, 61, 55, 49, 49, 59, 50, NA, 58, 45, NA, 55, 52, 64, NA, 49, 58, 57, 52, 60, NA, 61, NA, 31, 58, 52, 56, 5…##Part 1: Running our “Final” 2-level Model…DV=gktmathss (math scores), Student level IV: gkselfconcraw (self-concept), gkclasstype (classroom type) and gkthighdegree (teacher highest degree) as teacher/classroom level IVs.
model.final <- lmer(gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1|gktchid), REML = FALSE, data = projectSTAR)
summary(model.final)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1 | gktchid)
Data: projectSTAR
AIC BIC logLik deviance df.resid
49055.9 49094.8 -24522.0 49043.9 4750
Scaled residuals:
Min 1Q Median 3Q Max
-3.8391 -0.6386 -0.0618 0.5814 3.8970
Random effects:
Groups Name Variance Std.Dev.
gktchid (Intercept) 663.4 25.76
Residual 1550.2 39.37
Number of obs: 4756, groups: gktchid, 301
Fixed effects:
Estimate Std. Error t value
(Intercept) 464.8157 10.1601 45.749
gkselfconcraw 0.4458 0.1196 3.726
gkclasstype -4.6260 1.9336 -2.392
gkthighdegree 2.7561 2.7143 1.015
Correlation of Fixed Effects:
(Intr) gkslfc gkclss
gkselfcncrw -0.660
gkclasstype -0.377 0.017
gkthighdegr -0.630 -0.010 -0.002##Part 2: Begin to check assumptions ##Question 2: Use the “augment” function from the broom.mixed Package to Get Predictions, Residuals, Cook’s Distances, Etc.:
library(broom.mixed)
diagnostics <- augment(model.final)##Part 2 ##Question 2 cont’d: Assess Normality of Residuals (Visually, with a Histogram)
ggplot(data = diagnostics, mapping = aes(x = .resid)) +
geom_histogram(binwidth = .25) + theme_classic() +
labs(title = "Histogram of Residuals for Math Scores Model",
x = "Residual Value") +
geom_vline(xintercept = c(-2.5, 2.5), linetype = "dotted")
##Part 2 ##Question 2 cont’d: Shapiro Wilk Test
shapiro.test(diagnostics$.resid)
Shapiro-Wilk normality test
data: diagnostics$.resid
W = 0.99119, p-value < 2.2e-16##Question 3: Use Residuals vs. Fitted (RVF) Plot to Assess Homoskedasticity of Errors (Basic RVF Plot, Without Groups)
ggplot(data = diagnostics, mapping = aes(x = .fitted, y = .resid)) +
geom_point() + labs(title = "RVF Plot for Math Scores Model",
x = "Predicted Value, Math Scores",
y = "Residual Value") + theme_classic()
##trying same thing with color…???
ggplot(data = diagnostics, mapping = aes(x = .fitted, y = .resid, color = gktchid)) +
geom_point() + labs(title = "RVF Plot for Math Scores Model, by Kindergarten ID",
x = "Predicted Value, Math Scores",
y = "Residual Value") + theme_classic()
##attempt to make gktchid more meaningful…(is that even possible?)
math.clean <- projectSTAR %>%
mutate(.,
gktchid_1000 = gktchid/10000)
glimpse(math.clean)Rows: 6,325
Columns: 28
$ stdntid <dbl> 10001, 10133, 10246, 10263, 10266, 10275, 10281, 10282, 10285, 10286, 10287, 10289, 10290, 10291, 10292, 10294, 10295, 10299, 10300, 10302, 1…
$ race <dbl+lbl> 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2,…
$ gender <dbl+lbl> 1, 1, 2, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 1,…
$ FLAGSGK <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ flaggk <dbl+lbl> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,…
$ gkclasstype <dbl+lbl> 3, 3, 3, 1, 2, 3, 1, 3, 1, 2, 3, 3, 2, 3, 3, 3, 1, 2, 3, 3, 2, 1, 1, 3, 2, 1, 1, 2, 3, 2, 3, 3, 3, 2, 3, 2, 3, 3, 1, 1, 3, 1, 3, 1, 3, 1,…
$ gkschid <dbl> 169229, 169280, 218562, 205492, 257899, 161176, 189382, 189382, 201449, 230612, 128068, 262937, 169219, 169229, 205491, 159171, 161176, 16821…
$ gksurban <dbl+lbl> 2, 2, 4, 2, 3, 3, 2, 2, 3, 3, 3, 2, 2, 2, 2, 3, 3, 3, 3, 3, 2, 1, 2, 3, 1, 4, 2, 3, 2, 3, 4, 3, 2, 2, 2, 3, 2, 2, 4, 3, 1, 1, 1, 1, 1, 1,…
$ gktchid <dbl> 16922904, 16928003, 21856202, 20549204, 25789904, 16117602, 18938204, 18938203, 20144901, 23061203, 12806803, 26293701, 16921902, 16922904, 2…
$ gktgen <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
$ gktrace <dbl+lbl> 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, NA, 1, 1, 1, …
$ gkthighdegree <dbl+lbl> 2, 2, 3, 2, 3, 2, 3, 3, 3, 3, 2, 2, 3, 2, 2, 2, 2, 2, 2, 3, 2, 2, 3, 2, 2, 3, 3, 2, 2, 2, NA, 2, 2, 3, …
$ gktcareer <dbl+lbl> 4, 4, 4, 4, 4, NA, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, NA, 2, 4, 4, 2, 4, 4, 4, 4, NA, NA, 4, 4, 6, …
$ gktyears <dbl> 5, 7, 8, 3, 12, 2, 7, 14, 4, 6, 11, 16, 12, 5, 17, 10, 6, 10, 13, 9, 18, 1, 7, 13, 1, 14, 18, 8, 8, 4, NA, 13, 17, 15, 16, 13, 3, 17, 11, 2, …
$ gkclasssize <dbl> 24, 22, 17, 17, 24, 24, 13, 24, 14, 23, 23, 22, 20, 24, 23, 27, 17, 24, 22, 23, 22, 17, 13, 18, 20, 14, 14, 23, 24, 24, 21, 18, 23, 19, 22, 2…
$ gkfreelunch <dbl+lbl> 2, 2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 2, 1, 2, 2, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1,…
$ gkrepeat <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ gkspeced <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 1, 2, 2,…
$ gkspecin <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
$ gkpresent <dbl> 161, 175, 132, 178, 170, 94, 160, 154, 172, 95, 163, 172, 180, 149, 173, NA, 160, 173, 170, 161, 128, 151, 64, 84, 145, 160, 161, 177, 154, 1…
$ gkabsent <dbl> 19, 5, 28, 2, 10, 3, 2, 7, 8, 2, 17, 8, 0, 31, 7, NA, 20, 7, 2, 19, 24, 0, 8, 9, 7, 4, 19, 3, 3, 4, 15, 26, 6, 5, 10, 3, 15, 13, 2, 1, 12, 7,…
$ gktreadss <dbl> NA, 427, 450, 483, 456, 411, 443, 448, 463, 472, 428, 545, 408, 422, 472, NA, 437, 450, 497, 460, 408, NA, 467, 403, 411, 490, 503, NA, 442, …
$ gktmathss <dbl> NA, 478, 494, 513, 513, 468, 473, 449, 520, 536, 484, 626, 454, 439, 528, NA, 463, 484, 547, 576, 449, NA, 576, 392, 444, 559, 576, NA, 547, …
$ gktlistss <dbl> NA, 509, 549, 554, 520, 571, 595, 540, 565, 595, NA, 622, 474, 536, 578, NA, 554, 565, 524, 595, 512, NA, 565, 501, 505, 595, 671, NA, 532, 5…
$ gkwordskillss <dbl> NA, 418, 444, 431, 468, 396, 444, 444, 480, 486, 423, 524, 410, 423, 458, NA, 427, 435, 512, 449, 405, NA, 480, 410, 405, 512, 480, NA, 435, …
$ gkmotivraw <dbl> 23, 24, 28, 27, 25, 24, NA, NA, 26, 27, 24, 24, 23, 28, 24, NA, 26, 25, NA, 26, 25, 29, NA, 24, 26, 24, 25, 28, NA, 27, NA, 17, 24, 23, 25, 2…
$ gkselfconcraw <dbl> 52, 53, 56, 61, 54, 55, NA, NA, 52, 61, 55, 49, 49, 59, 50, NA, 58, 45, NA, 55, 52, 64, NA, 49, 58, 57, 52, 60, NA, 61, NA, 31, 58, 52, 56, 5…
$ gktchid_1000 <dbl> 1692.290, 1692.800, 2185.620, 2054.920, 2578.990, 1611.760, 1893.820, 1893.820, 2014.490, 2306.120, 1280.680, 2629.370, 1692.190, 1692.290, 2…##redo the color plot to see if it is easier to see what is going…
ggplot(data = diagnostics, mapping = aes(x = .fitted, y = .resid, color = gktchid_1000)) +
geom_point() + labs(title = "RVF Plot for Math Scores Model, by Kindergarten ID",
x = "Predicted Value, Math Scores",
y = "Residual Value") + theme_classic()##Oh I don’t think it is working because I did not include this divided variable in my original “final model” Moving on :)
##Question 4: Using Cook’s Distance to Identify Outliers
ggplot(data = diagnostics, mapping = aes(x = .fitted, y = .cooksd, label = gktchid)) +
geom_point() + geom_text(nudge_x = .25) + theme_classic() +
labs(title = "Cook's Distance Plot for Math Scores Model",
x = "Predicted Value, % Math Scores",
y = "Cook's Distance") +
geom_hline(yintercept = 4/816, linetype = "dotted")
##Part 3 ##Robustness check! AKA:re-Run Model with Robust Standard Errors
library(robustlmm)
model.robust <- rlmer(gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1|gktchid), REML=FALSE, data = projectSTAR)
summary(model.robust)Robust linear mixed model fit by DAStau
Formula: gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1 | gktchid)
Data: projectSTAR
Scaled residuals:
Min 1Q Median 3Q Max
-4.1521 -0.6275 -0.0298 0.6357 4.1638
Random effects:
Groups Name Variance Std.Dev.
gktchid (Intercept) 619.8 24.90
Residual 1396.2 37.37
Number of obs: 4756, groups: gktchid, 301
Fixed effects:
Estimate Std. Error t value
(Intercept) 461.6699 9.9800 46.26
gkselfconcraw 0.4537 0.1165 3.89
gkclasstype -4.0642 1.9121 -2.13
gkthighdegree 2.7485 2.6839 1.02
Correlation of Fixed Effects:
(Intr) gkslfc gkclss
gkselfcncrw -0.654
gkclasstype -0.379 0.017
gkthighdegr -0.635 -0.009 -0.002
Robustness weights for the residuals:
3806 weights are ~= 1. The remaining 950 ones are summarized as
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.323 0.642 0.809 0.775 0.937 0.999
Robustness weights for the random effects:
241 weights are ~= 1. The remaining 60 ones are summarized as
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.353 0.667 0.889 0.810 0.971 0.998
Rho functions used for fitting:
Residuals:
eff: smoothed Huber (k = 1.345, s = 10)
sig: smoothed Huber, Proposal II (k = 1.345, s = 10)
Random Effects, variance component 1 (gktchid):
eff: smoothed Huber (k = 1.345, s = 10)
vcp: smoothed Huber, Proposal II (k = 1.345, s = 10) ##re-run model with outliers removed (Cook’s Distance >.10)
STAR.trimmed <- diagnostics %>%
filter(., .cooksd < .10)
model.trimmed <- lmer(gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1|gktchid), REML=FALSE, data = STAR.trimmed)
summary(model.trimmed)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1 | gktchid)
Data: STAR.trimmed
AIC BIC logLik deviance df.resid
47346.0 47384.7 -23667.0 47334.0 4657
Scaled residuals:
Min 1Q Median 3Q Max
-3.05922 -0.65534 -0.04224 0.61962 3.08934
Random effects:
Groups Name Variance Std.Dev.
gktchid (Intercept) 704.3 26.54
Residual 1300.7 36.06
Number of obs: 4663, groups: gktchid, 301
Fixed effects:
Estimate Std. Error t value
(Intercept) 461.5853 9.9617 46.336
gkselfconcraw 0.4567 0.1113 4.102
gkclasstype -4.0598 1.9670 -2.064
gkthighdegree 2.8253 2.7594 1.024
Correlation of Fixed Effects:
(Intr) gkslfc gkclss
gkselfcncrw -0.627
gkclasstype -0.390 0.017
gkthighdegr -0.655 -0.009 -0.002##Now time to make a table to compare original “final” model, Robust Model and Trimmed Model for Differences (sensetivity analysis)
library(modelsummary)
models <- list(model.final, model.trimmed)
modelsummary(models, output = "markdown")| Model 1 | Model 2 | |
|---|---|---|
| (Intercept) | 464.816 | 461.585 |
| (10.160) | (9.962) | |
| gkselfconcraw | 0.446 | 0.457 |
| (0.120) | (0.111) | |
| gkclasstype | -4.626 | -4.060 |
| (1.934) | (1.967) | |
| gkthighdegree | 2.756 | 2.825 |
| (2.714) | (2.759) | |
| sd__(Intercept) | 25.756 | 26.539 |
| sd__Observation | 39.373 | 36.065 |
| AIC | 49055.9 | 47346.0 |
| BIC | 49094.8 | 47384.7 |
| Log.Lik. | -24521.975 | -23667.004 |
LS0tCnRpdGxlOiAiTW9kdWxlIDEwIgphdXRob3I6ICJDYWl0bGluIE1heXRvbiIKb3V0cHV0OiBodG1sX25vdGVib29rCgotLS0KClRoaXMgd2VlayB3ZSB1c2VkIFByb2plY3RTVEFSIGRhdGEgYWdhaW4sIHdoaWNoIGlzLi4uImEgc2FtcGxlIG9mIG1vcmUgdGhhbiA2LDAwMCBzdHVkZW50cyBuZXN0ZWQgd2l0aGluIDIwMCsgc2Nob29scyBjb2xsZWN0ZWQgYmV0d2VlbiAxOTg1IGFuZCAxOTkwIGFzIHBhcnQgb2YgdGhlIFRlbm5lc3NlZSBDbGFzcyBTaXplIEV4cGVyaW1lbnQuIiAtRHIuIEJyb2RhIChNb2R1bGUgNykKICAKCiNMb2FkIGluIG91ciBNVlAgcGFja2FnZXMKCmBgYHtyfQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMobGlicmFyeSh0aWR5dmVyc2UpKQpsaWJyYXJ5KGxtZTQpCmxpYnJhcnkoSG1pc2MpCmBgYAoKIyBQYXJ0IDE6IGxvYWRpbmcgaW4gdGhlIGRhdGEKYGBge3J9Cgpwcm9qZWN0U1RBUiA8LSBoYXZlbjo6cmVhZF9kdGEoInByb2plY3RTVEFSLmR0YSIpCgpnbGltcHNlKHByb2plY3RTVEFSKQoKYGBgCgojI1BhcnQgMTogUnVubmluZyBvdXIgIkZpbmFsIiAyLWxldmVsIE1vZGVsLi4uRFY9Z2t0bWF0aHNzIChtYXRoIHNjb3JlcyksIFN0dWRlbnQgbGV2ZWwgSVY6IGdrc2VsZmNvbmNyYXcgKHNlbGYtY29uY2VwdCksIGdrY2xhc3N0eXBlIChjbGFzc3Jvb20gdHlwZSkgYW5kIGdrdGhpZ2hkZWdyZWUgKHRlYWNoZXIgaGlnaGVzdCBkZWdyZWUpIGFzIHRlYWNoZXIvY2xhc3Nyb29tIGxldmVsIElWcy4gCmBgYHtyfQptb2RlbC5maW5hbCA8LSBsbWVyKGdrdG1hdGhzcyB+IGdrc2VsZmNvbmNyYXcgKyBna2NsYXNzdHlwZSArIGdrdGhpZ2hkZWdyZWUgKyAoMXxna3RjaGlkKSwgUkVNTCA9IEZBTFNFLCBkYXRhID0gcHJvamVjdFNUQVIpCnN1bW1hcnkobW9kZWwuZmluYWwpCmBgYAoKIyNQYXJ0IDI6IEJlZ2luIHRvIGNoZWNrIGFzc3VtcHRpb25zCiMjUXVlc3Rpb24gMjogVXNlIHRoZSAiYXVnbWVudCIgZnVuY3Rpb24gZnJvbSB0aGUgYnJvb20ubWl4ZWQgUGFja2FnZSB0byBHZXQgUHJlZGljdGlvbnMsIFJlc2lkdWFscywgQ29va+KAmXMgRGlzdGFuY2VzLCBFdGMuOgpgYGB7cn0KbGlicmFyeShicm9vbS5taXhlZCkKCmRpYWdub3N0aWNzIDwtIGF1Z21lbnQobW9kZWwuZmluYWwpCmBgYAoKIyNQYXJ0IDIKIyNRdWVzdGlvbiAyIGNvbnQnZDogQXNzZXNzIE5vcm1hbGl0eSBvZiBSZXNpZHVhbHMgKFZpc3VhbGx5LCB3aXRoIGEgSGlzdG9ncmFtKQpgYGB7cn0KZ2dwbG90KGRhdGEgPSBkaWFnbm9zdGljcywgbWFwcGluZyA9IGFlcyh4ID0gLnJlc2lkKSkgKwogIGdlb21faGlzdG9ncmFtKGJpbndpZHRoID0gLjI1KSArIHRoZW1lX2NsYXNzaWMoKSArIAogIGxhYnModGl0bGUgPSAiSGlzdG9ncmFtIG9mIFJlc2lkdWFscyBmb3IgTWF0aCBTY29yZXMgTW9kZWwiLAogICAgICAgICAgICAgICAgICAgICAgeCA9ICJSZXNpZHVhbCBWYWx1ZSIpICsKICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSBjKC0yLjUsIDIuNSksIGxpbmV0eXBlID0gImRvdHRlZCIpCmBgYAojI1BhcnQgMgojI1F1ZXN0aW9uIDIgY29udCdkOiBTaGFwaXJvIFdpbGsgVGVzdApgYGB7cn0Kc2hhcGlyby50ZXN0KGRpYWdub3N0aWNzJC5yZXNpZCkKYGBgCgojI1F1ZXN0aW9uIDM6IFVzZSBSZXNpZHVhbHMgdnMuIEZpdHRlZCAoUlZGKSBQbG90IHRvIEFzc2VzcyBIb21vc2tlZGFzdGljaXR5IG9mIEVycm9ycyAoQmFzaWMgUlZGIFBsb3QsIFdpdGhvdXQgR3JvdXBzKQpgYGB7cn0KZ2dwbG90KGRhdGEgPSBkaWFnbm9zdGljcywgbWFwcGluZyA9IGFlcyh4ID0gLmZpdHRlZCwgeSA9IC5yZXNpZCkpICsKICBnZW9tX3BvaW50KCkgKyBsYWJzKHRpdGxlID0gIlJWRiBQbG90IGZvciBNYXRoIFNjb3JlcyBNb2RlbCIsCiAgICAgICAgICAgICAgICAgICAgICB4ID0gIlByZWRpY3RlZCBWYWx1ZSwgTWF0aCBTY29yZXMiLAogICAgICAgICAgICAgICAgICAgICAgeSA9ICJSZXNpZHVhbCBWYWx1ZSIpICsgdGhlbWVfY2xhc3NpYygpCmBgYAoKIyN0cnlpbmcgc2FtZSB0aGluZyB3aXRoIGNvbG9yLi4uPz8/CmBgYHtyfQpnZ3Bsb3QoZGF0YSA9IGRpYWdub3N0aWNzLCBtYXBwaW5nID0gYWVzKHggPSAuZml0dGVkLCB5ID0gLnJlc2lkLCBjb2xvciA9IGdrdGNoaWQpKSArCiAgZ2VvbV9wb2ludCgpICsgbGFicyh0aXRsZSA9ICJSVkYgUGxvdCBmb3IgTWF0aCBTY29yZXMgTW9kZWwsIGJ5IEtpbmRlcmdhcnRlbiBJRCIsCiAgICAgICAgICAgICAgICAgICAgICB4ID0gIlByZWRpY3RlZCBWYWx1ZSwgTWF0aCBTY29yZXMiLAogICAgICAgICAgICAgICAgICAgICAgeSA9ICJSZXNpZHVhbCBWYWx1ZSIpICsgdGhlbWVfY2xhc3NpYygpCmBgYAojI2F0dGVtcHQgdG8gbWFrZSBna3RjaGlkIG1vcmUgbWVhbmluZ2Z1bC4uLihpcyB0aGF0IGV2ZW4gcG9zc2libGU/KQoKYGBge3J9Cm1hdGguY2xlYW4gPC0gcHJvamVjdFNUQVIgJT4lCiAgbXV0YXRlKC4sCiAgICAgICAgIGdrdGNoaWRfMTAwMCA9IGdrdGNoaWQvMTAwMDApCgpnbGltcHNlKG1hdGguY2xlYW4pCmBgYAoKIyNyZWRvIHRoZSBjb2xvciBwbG90IHRvIHNlZSBpZiBpdCBpcyBlYXNpZXIgdG8gc2VlIHdoYXQgaXMgZ29pbmcuLi4KYGBge3J9CmdncGxvdChkYXRhID0gZGlhZ25vc3RpY3MsIG1hcHBpbmcgPSBhZXMoeCA9IC5maXR0ZWQsIHkgPSAucmVzaWQsIGNvbG9yID0gZ2t0Y2hpZF8xMDAwKSkgKwogIGdlb21fcG9pbnQoKSArIGxhYnModGl0bGUgPSAiUlZGIFBsb3QgZm9yIE1hdGggU2NvcmVzIE1vZGVsLCBieSBLaW5kZXJnYXJ0ZW4gSUQiLAogICAgICAgICAgICAgICAgICAgICAgeCA9ICJQcmVkaWN0ZWQgVmFsdWUsIE1hdGggU2NvcmVzIiwKICAgICAgICAgICAgICAgICAgICAgIHkgPSAiUmVzaWR1YWwgVmFsdWUiKSArIHRoZW1lX2NsYXNzaWMoKQpgYGAKCiMjT2ggSSBkb24ndCB0aGluayBpdCBpcyB3b3JraW5nIGJlY2F1c2UgSSBkaWQgbm90IGluY2x1ZGUgdGhpcyBkaXZpZGVkIHZhcmlhYmxlIGluIG15IG9yaWdpbmFsICJmaW5hbCBtb2RlbCIgTW92aW5nIG9uIDopCgojI1F1ZXN0aW9uIDQ6IFVzaW5nIENvb2vigJlzIERpc3RhbmNlIHRvIElkZW50aWZ5IE91dGxpZXJzCmBgYHtyfQpnZ3Bsb3QoZGF0YSA9IGRpYWdub3N0aWNzLCBtYXBwaW5nID0gYWVzKHggPSAuZml0dGVkLCB5ID0gLmNvb2tzZCwgbGFiZWwgPSBna3RjaGlkKSkgKwogIGdlb21fcG9pbnQoKSArIGdlb21fdGV4dChudWRnZV94ID0gLjI1KSArIHRoZW1lX2NsYXNzaWMoKSArIAogIGxhYnModGl0bGUgPSAiQ29vaydzIERpc3RhbmNlIFBsb3QgZm9yIE1hdGggU2NvcmVzIE1vZGVsIiwKICAgICAgICAgICAgICAgICAgICAgIHggPSAiUHJlZGljdGVkIFZhbHVlLCAlIE1hdGggU2NvcmVzIiwKICAgICAgICAgICAgICAgICAgICAgIHkgPSAiQ29vaydzIERpc3RhbmNlIikgKyAKICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSA0LzgxNiwgbGluZXR5cGUgPSAiZG90dGVkIikKYGBgCgojI1BhcnQgMwojI1JvYnVzdG5lc3MgY2hlY2shIEFLQTpyZS1SdW4gTW9kZWwgd2l0aCBSb2J1c3QgU3RhbmRhcmQgRXJyb3JzCmBgYHtyfQpsaWJyYXJ5KHJvYnVzdGxtbSkKbW9kZWwucm9idXN0IDwtIHJsbWVyKGdrdG1hdGhzcyB+IGdrc2VsZmNvbmNyYXcgKyBna2NsYXNzdHlwZSArIGdrdGhpZ2hkZWdyZWUgKyAoMXxna3RjaGlkKSwgUkVNTD1GQUxTRSwgZGF0YSA9IHByb2plY3RTVEFSKQpzdW1tYXJ5KG1vZGVsLnJvYnVzdCkKYGBgCgojI3JlLXJ1biBtb2RlbCB3aXRoIG91dGxpZXJzIHJlbW92ZWQgKENvb2sncyBEaXN0YW5jZSA+LjEwKQpgYGB7cn0KU1RBUi50cmltbWVkIDwtIGRpYWdub3N0aWNzICU+JQogIGZpbHRlciguLCAuY29va3NkIDwgLjEwKQptb2RlbC50cmltbWVkIDwtIGxtZXIoZ2t0bWF0aHNzIH4gZ2tzZWxmY29uY3JhdyArIGdrY2xhc3N0eXBlICsgZ2t0aGlnaGRlZ3JlZSArICgxfGdrdGNoaWQpLCBSRU1MPUZBTFNFLCBkYXRhID0gU1RBUi50cmltbWVkKQpzdW1tYXJ5KG1vZGVsLnRyaW1tZWQpCmBgYAoKIyNOb3cgdGltZSB0byBtYWtlIGEgdGFibGUgdG8gY29tcGFyZSBvcmlnaW5hbCAiZmluYWwiIG1vZGVsLCBSb2J1c3QgTW9kZWwgYW5kIFRyaW1tZWQgTW9kZWwgZm9yIERpZmZlcmVuY2VzIChzZW5zZXRpdml0eSBhbmFseXNpcykKYGBge3J9CmxpYnJhcnkobW9kZWxzdW1tYXJ5KQptb2RlbHMgPC0gbGlzdChtb2RlbC5maW5hbCwgbW9kZWwudHJpbW1lZCkKbW9kZWxzdW1tYXJ5KG1vZGVscywgb3V0cHV0ID0gIm1hcmtkb3duIikKYGBgCgo=