Multilevel Modeling, Module 10: Checking Assumptions for MLMs
Dr. Broda
This week, we will us the productivity dataset from Stata, developed by Dr. Alicia Munnell and published by the Federal Reserve Bank of Boston (1990).
Load in Our MVP Packages
suppressPackageStartupMessages(library(tidyverse))
library(lme4)
library(Hmisc)Part 1: Loading in the Data, and then Recoding Variables
Load in the Data
productivity <- haven::read_dta("productivity.dta")
glimpse(productivity)Rows: 816
Columns: 11
$ state [3m[38;5;246m<dbl>[39m[23m 19, 44, 5, 41, 45, 16, 37, 20, 25, 11, 14, 40, 23, 42, 4, 26, 47, 2, 33, 43, 28, 6, 39, 7, 30, 38, …
$ region [3m[38;5;246m<dbl>[39m[23m 1, 5, 8, 7, 9, 7, 1, 3, 4, 3, 4, 6, 4, 8, 9, 8, 3, 8, 3, 1, 2, 1, 4, 5, 2, 5, 5, 7, 4, 6, 1, 5, 6, …
$ year [3m[38;5;246m<dbl>[39m[23m 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 197…
$ public [3m[38;5;246m<dbl>[39m[23m 21923.28, 20732.56, 11402.52, 53639.88, 25751.49, 19638.31, 4310.22, 44684.82, 9235.03, 52197.49, 1…
$ hwy [3m[38;5;246m<dbl>[39m[23m 9.052737, 9.329386, 8.390089, 10.163615, 9.066138, 9.199419, 7.725149, 9.793947, 8.327134, 10.06827…
$ water [3m[38;5;246m<dbl>[39m[23m 7.714320, 7.892493, 7.680190, 9.030341, 8.142238, 7.886390, 6.349594, 9.033661, 6.971622, 8.758089,…
$ other [3m[38;5;246m<dbl>[39m[23m 9.318403, 8.823389, 8.483487, 9.870202, 9.522063, 8.865893, 7.295308, 9.818987, 8.302643, 10.010218…
$ private [3m[38;5;246m<dbl>[39m[23m 10.543170, 10.454692, 10.073642, 12.272218, 10.529598, 11.490705, 8.601034, 11.347531, 9.860606, 11…
$ gsp [3m[38;5;246m<dbl>[39m[23m 11.073319, 10.784545, 10.153818, 11.988675, 10.577044, 11.045574, 9.192584, 11.534413, 9.684523, 11…
$ emp [3m[38;5;246m<dbl>[39m[23m 7.715793, 7.325742, 6.620340, 8.195582, 6.984160, 6.940803, 5.840932, 8.006034, 6.182704, 8.376919,…
$ unemp [3m[38;5;246m<dbl>[39m[23m 4.6, 3.4, 4.4, 4.4, 9.1, 6.6, 5.2, 6.7, 3.1, 4.1, 4.8, 4.8, 3.3, 6.1, 7.2, 5.9, 3.9, 4.4, 5.4, 4.9,…The describe function from Hmisc is a nice codebook…
Hmisc::describe(productivity)productivity
11 Variables 816 Observations
----------------------------------------------------------------------------------------------------------------------
state : states 1-48 Format:%9.0g
n missing distinct Info Mean Gmd .05 .10 .25 .50 .75 .90 .95
816 0 48 1 24.5 16.01 3.00 5.00 12.75 24.50 36.25 44.00 46.00
lowest : 1 2 3 4 5, highest: 44 45 46 47 48
----------------------------------------------------------------------------------------------------------------------
region : regions 1-9 Format:%9.0g
n missing distinct Info Mean Gmd
816 0 9 0.983 4.958 2.811
lowest : 1 2 3 4 5, highest: 5 6 7 8 9
Value 1 2 3 4 5 6 7 8 9
Frequency 102 51 85 119 136 68 68 136 51
Proportion 0.125 0.062 0.104 0.146 0.167 0.083 0.083 0.167 0.062
----------------------------------------------------------------------------------------------------------------------
year : years 1970-1986 Format:%9.0g
n missing distinct Info Mean Gmd .05 .10 .25 .50 .75 .90 .95
816 0 17 0.997 1978 5.654 1970 1971 1974 1978 1982 1985 1986
lowest : 1970 1971 1972 1973 1974, highest: 1982 1983 1984 1985 1986
Value 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986
Frequency 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48
Proportion 0.059 0.059 0.059 0.059 0.059 0.059 0.059 0.059 0.059 0.059 0.059 0.059 0.059 0.059 0.059 0.059 0.059
----------------------------------------------------------------------------------------------------------------------
public : public capital stock Format:%9.0g
n missing distinct Info Mean Gmd .05 .10 .25 .50 .75 .90 .95
816 0 815 1 25037 24904 4025 4407 7097 17572 27692 56488 69628
lowest : 2627.12 2756.19 2894.58 2925.39 2927.85, highest: 139514.19 139662.44 139776.97 139910.22 140217.31
----------------------------------------------------------------------------------------------------------------------
hwy : log(highway component of public) Format:%9.0g
n missing distinct Info Mean Gmd .05 .10 .25 .50 .75 .90 .95
816 0 816 1 8.9 0.9192 7.639 7.758 8.258 8.930 9.330 10.113 10.308
lowest : 7.510507 7.512159 7.516809 7.521551 7.526922, highest: 10.759274 10.760835 10.769929 10.770011 10.772675
----------------------------------------------------------------------------------------------------------------------
water : log(water component of public) Format:%9.0g
n missing distinct Info Mean Gmd .05 .10 .25 .50 .75 .90 .95
816 0 815 1 7.6 1.3 5.662 5.962 6.639 7.726 8.371 9.097 9.299
lowest : 5.431361 5.437470 5.440294 5.453568 5.456304, highest: 10.064859 10.078645 10.079176 10.090852 10.110189
----------------------------------------------------------------------------------------------------------------------
other : log(bldg/other component of public) Format:%9.0g
n missing distinct Info Mean Gmd .05 .10 .25 .50 .75 .90 .95
816 0 816 1 8.736 1.245 6.995 7.202 7.819 8.855 9.359 10.067 10.265
lowest : 6.288769 6.427103 6.505784 6.507755 6.515690, highest: 11.275786 11.289294 11.296843 11.297787 11.298842
----------------------------------------------------------------------------------------------------------------------
private : log(private capital stock) Format:%9.0g
n missing distinct Info Mean Gmd .05 .10 .25 .50 .75 .90 .95
816 0 816 1 10.56 1.044 8.926 9.285 9.983 10.613 11.079 11.791 12.061
lowest : 8.307141 8.354745 8.390010 8.428930 8.474315, highest: 12.791270 12.799979 12.804304 12.805800 12.835592
----------------------------------------------------------------------------------------------------------------------
gsp : log(gross state product) Format:%9.0g
n missing distinct Info Mean Gmd .05 .10 .25 .50 .75 .90 .95
816 0 807 1 10.51 1.165 8.961 9.099 9.711 10.596 11.129 11.882 12.392
lowest : 8.378850 8.397959 8.418036 8.436200 8.445052, highest: 12.848508 12.875255 12.949259 13.003764 13.048824
----------------------------------------------------------------------------------------------------------------------
emp : log(non-agriculture payrolls) Format:%9.0g
n missing distinct Info Mean Gmd .05 .10 .25 .50 .75 .90 .95
816 0 807 1 6.978 1.166 5.332 5.516 6.163 7.060 7.656 8.381 8.642
lowest : 4.684905 4.709530 4.764735 4.837075 4.916325, highest: 9.206915 9.208869 9.266153 9.303740 9.328835
----------------------------------------------------------------------------------------------------------------------
unemp : state unemployment rate Format:%9.0g
n missing distinct Info Mean Gmd .05 .10 .25 .50 .75 .90 .95
816 0 80 1 6.602 2.456 3.6 4.0 5.0 6.2 7.9 9.5 11.0
lowest : 2.8 2.9 3.0 3.1 3.2, highest: 13.0 14.0 15.0 16.0 18.0
----------------------------------------------------------------------------------------------------------------------A little quick cleaning - Convert Region to Factor Variable, Rescale Public Capital Stock, “Zero Out” Year and Create Nonlinear Terms
prod.clean <- productivity %>%
mutate(.,
region.fac = as_factor(region),
year0 = year - 1970,
year_sq = year0^2,
year_cube = year0^3,
pub_cap_1000 = public/1000)
glimpse(prod.clean)Rows: 816
Columns: 16
$ state [3m[38;5;246m<dbl>[39m[23m 19, 44, 5, 41, 45, 16, 37, 20, 25, 11, 14, 40, 23, 42, 4, 26, 47, 2, 33, 43, 28, 6, 39, 7, 30,…
$ region [3m[38;5;246m<dbl>[39m[23m 1, 5, 8, 7, 9, 7, 1, 3, 4, 3, 4, 6, 4, 8, 9, 8, 3, 8, 3, 1, 2, 1, 4, 5, 2, 5, 5, 7, 4, 6, 1, 5…
$ year [3m[38;5;246m<dbl>[39m[23m 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970, 1970…
$ public [3m[38;5;246m<dbl>[39m[23m 21923.28, 20732.56, 11402.52, 53639.88, 25751.49, 19638.31, 4310.22, 44684.82, 9235.03, 52197.…
$ hwy [3m[38;5;246m<dbl>[39m[23m 9.052737, 9.329386, 8.390089, 10.163615, 9.066138, 9.199419, 7.725149, 9.793947, 8.327134, 10.…
$ water [3m[38;5;246m<dbl>[39m[23m 7.714320, 7.892493, 7.680190, 9.030341, 8.142238, 7.886390, 6.349594, 9.033661, 6.971622, 8.75…
$ other [3m[38;5;246m<dbl>[39m[23m 9.318403, 8.823389, 8.483487, 9.870202, 9.522063, 8.865893, 7.295308, 9.818987, 8.302643, 10.0…
$ private [3m[38;5;246m<dbl>[39m[23m 10.543170, 10.454692, 10.073642, 12.272218, 10.529598, 11.490705, 8.601034, 11.347531, 9.86060…
$ gsp [3m[38;5;246m<dbl>[39m[23m 11.073319, 10.784545, 10.153818, 11.988675, 10.577044, 11.045574, 9.192584, 11.534413, 9.68452…
$ emp [3m[38;5;246m<dbl>[39m[23m 7.715793, 7.325742, 6.620340, 8.195582, 6.984160, 6.940803, 5.840932, 8.006034, 6.182704, 8.37…
$ unemp [3m[38;5;246m<dbl>[39m[23m 4.6, 3.4, 4.4, 4.4, 9.1, 6.6, 5.2, 6.7, 3.1, 4.1, 4.8, 4.8, 3.3, 6.1, 7.2, 5.9, 3.9, 4.4, 5.4,…
$ region.fac [3m[38;5;246m<fct>[39m[23m 1, 5, 8, 7, 9, 7, 1, 3, 4, 3, 4, 6, 4, 8, 9, 8, 3, 8, 3, 1, 2, 1, 4, 5, 2, 5, 5, 7, 4, 6, 1, 5…
$ year0 [3m[38;5;246m<dbl>[39m[23m 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ year_sq [3m[38;5;246m<dbl>[39m[23m 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ year_cube [3m[38;5;246m<dbl>[39m[23m 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ pub_cap_1000 [3m[38;5;246m<dbl>[39m[23m 21.92328, 20.73256, 11.40252, 53.63988, 25.75149, 19.63831, 4.31022, 44.68482, 9.23503, 52.197…Part 2: Find Your “Final” or Preferred Model
Start with a null growth model since this is longitudinal data…
model.null.growth <- lmer(unemp ~ year0 + (1|state), REML = FALSE, data = prod.clean)
summary(model.null.growth)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: unemp ~ year0 + (1 | state)
Data: prod.clean
AIC BIC logLik deviance df.resid
3268.6 3287.4 -1630.3 3260.6 812
Scaled residuals:
Min 1Q Median 3Q Max
-2.6688 -0.6260 -0.1167 0.4891 4.8594
Random effects:
Groups Name Variance Std.Dev.
state (Intercept) 1.428 1.195
Residual 2.785 1.669
Number of obs: 816, groups: state, 48
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.17075 0.20557 25.15
year0 0.17893 0.01193 15.01
Correlation of Fixed Effects:
(Intr)
year0 -0.464Add in Quad Effects
model.2 <- lmer(unemp ~ year0 + year_sq + (1|state), REML = FALSE, data = prod.clean)
summary(model.2)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: unemp ~ year0 + year_sq + (1 | state)
Data: prod.clean
AIC BIC logLik deviance df.resid
3257.7 3281.2 -1623.8 3247.7 811
Scaled residuals:
Min 1Q Median 3Q Max
-2.4550 -0.6449 -0.1213 0.4901 4.9036
Random effects:
Groups Name Variance Std.Dev.
state (Intercept) 1.431 1.196
Residual 2.738 1.655
Number of obs: 816, groups: state, 48
Fixed effects:
Estimate Std. Error t value
(Intercept) 4.778638 0.232126 20.586
year0 0.335777 0.044987 7.464
year_sq -0.009803 0.002713 -3.614
Correlation of Fixed Effects:
(Intr) year0
year0 -0.558
year_sq 0.467 -0.965Add in Cubic Effects
model.3 <- lmer(unemp ~ year0 + year_cube + (1|state), REML = FALSE, data = prod.clean)Some predictor variables are on very different scales: consider rescalingsummary(model.3)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: unemp ~ year0 + year_cube + (1 | state)
Data: prod.clean
AIC BIC logLik deviance df.resid
3252.2 3275.7 -1621.1 3242.2 811
Scaled residuals:
Min 1Q Median 3Q Max
-2.3740 -0.6587 -0.1254 0.4910 4.8949
Random effects:
Groups Name Variance Std.Dev.
state (Intercept) 1.432 1.197
Residual 2.719 1.649
Number of obs: 816, groups: state, 48
Fixed effects:
Estimate Std. Error t value
(Intercept) 4.7913162 0.2230693 21.48
year0 0.2915168 0.0286038 10.19
year_cube -0.0004791 0.0001109 -4.32
Correlation of Fixed Effects:
(Intr) year0
year0 -0.533
year_cube 0.394 -0.911
fit warnings:
Some predictor variables are on very different scales: consider rescalingBoth Quadratic and Cubic Effects
model.4 <- lmer(unemp ~ year0 + year_sq + year_cube + (1|state), REML = FALSE, data = prod.clean)Some predictor variables are on very different scales: consider rescalingsummary(model.4)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: unemp ~ year0 + year_sq + year_cube + (1 | state)
Data: prod.clean
AIC BIC logLik deviance df.resid
3241.1 3269.3 -1614.5 3229.1 810
Scaled residuals:
Min 1Q Median 3Q Max
-2.2040 -0.6472 -0.1329 0.4740 4.8086
Random effects:
Groups Name Variance Std.Dev.
state (Intercept) 1.434 1.198
Residual 2.673 1.635
Number of obs: 816, groups: state, 48
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.2388480 0.2542657 20.604
year0 -0.0723854 0.1040941 -0.695
year_sq 0.0559415 0.0153967 3.633
year_cube -0.0027393 0.0006317 -4.336
Correlation of Fixed Effects:
(Intr) year0 yer_sq
year0 -0.592
year_sq 0.484 -0.962
year_cube -0.417 0.904 -0.985
fit warnings:
Some predictor variables are on very different scales: consider rescalingNow for some more predictors…
model.5 <- lmer(unemp ~ year0 + year_sq + year_cube + gsp + private + pub_cap_1000 + (1|state), REML = FALSE, data = prod.clean)Some predictor variables are on very different scales: consider rescalingsummary(model.5)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: unemp ~ year0 + year_sq + year_cube + gsp + private + pub_cap_1000 + (1 | state)
Data: prod.clean
AIC BIC logLik deviance df.resid
3108.9 3151.3 -1545.5 3090.9 807
Scaled residuals:
Min 1Q Median 3Q Max
-2.5959 -0.6424 -0.1207 0.5125 4.7506
Random effects:
Groups Name Variance Std.Dev.
state (Intercept) 9.703 3.115
Residual 1.993 1.412
Number of obs: 816, groups: state, 48
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.704e+01 5.335e+00 5.068
year0 -6.084e-03 9.176e-02 -0.066
year_sq 3.750e-02 1.338e-02 2.803
year_cube -1.787e-03 5.502e-04 -3.247
gsp -1.073e+01 7.365e-01 -14.574
private 8.409e+00 8.502e-01 9.891
pub_cap_1000 9.064e-02 1.686e-02 5.376
Correlation of Fixed Effects:
(Intr) year0 yer_sq yer_cb gsp privat
year0 0.134
year_sq 0.019 -0.944
year_cube -0.029 0.882 -0.984
gsp -0.127 -0.043 0.103 -0.126
private -0.518 -0.058 -0.089 0.116 -0.780
pub_cp_1000 0.550 -0.012 0.026 -0.011 -0.167 -0.232
fit warnings:
Some predictor variables are on very different scales: consider rescalingUse the augment function from the broom.mixed Package to Get Predictions, Residuals, Cook’s Distances, Etc.:
library(broom.mixed)Registered S3 methods overwritten by 'broom.mixed':
method from
augment.lme broom
augment.merMod broom
glance.lme broom
glance.merMod broom
glance.stanreg broom
tidy.brmsfit broom
tidy.gamlss broom
tidy.lme broom
tidy.merMod broom
tidy.rjags broom
tidy.stanfit broom
tidy.stanreg broomdiagnostics <- augment(model.5)
Step 1: 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 State Unemployment Model",
x = "Residual Value") +
geom_vline(xintercept = c(-2.5, 2.5), linetype = "dotted")
Statistically, with the Shapiro-Wilk Test
shapiro.test(diagnostics$.resid)
Shapiro-Wilk normality test
data: diagnostics$.resid
W = 0.97666, p-value = 3.899e-10Step 2: 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 State Unemployment Model",
x = "Predicted Value, % Unemployment",
y = "Residual Value") + theme_classic()
Same Thing, But With Added Color By Year0
ggplot(data = diagnostics, mapping = aes(x = .fitted, y = .resid, color = year0)) +
geom_point() + labs(title = "RVF Plot for State Unemployment Model, by Year",
x = "Predicted Value, % Unemployment",
y = "Residual Value") + theme_classic()
Step 3: Use Cook’s Distance to Identify Outliers
ggplot(data = diagnostics, mapping = aes(x = .fitted, y = .cooksd, label = state)) +
geom_point() + geom_text(nudge_x = .25) + theme_classic() +
labs(title = "Cook's Distance Plot for State Unemployment Model",
x = "Predicted Value, % Unemployment",
y = "Cook's Distance") +
geom_hline(yintercept = 4/816, linetype = "dotted")
NAStep 4: Re-Run Model with Robust Standard Errors
library(robustlmm)
model.robust <- rlmer(unemp ~ year0 + year_sq + year_cube + gsp + private + pub_cap_1000 + (1|state), REML = FALSE, data = prod.clean)Some predictor variables are on very different scales: consider rescalingsummary(model.robust)Robust linear mixed model fit by DAStau
Formula: unemp ~ year0 + year_sq + year_cube + gsp + private + pub_cap_1000 + (1 | state)
Data: prod.clean
Scaled residuals:
Min 1Q Median 3Q Max
-2.9180 -0.6434 -0.0609 0.5851 5.5981
Random effects:
Groups Name Variance Std.Dev.
state (Intercept) 7.072 2.659
Residual 1.729 1.315
Number of obs: 816, groups: state, 48
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.606e+01 4.862e+00 5.359
year0 4.230e-03 8.744e-02 0.048
year_sq 3.277e-02 1.278e-02 2.565
year_cube -1.512e-03 5.255e-04 -2.877
gsp -1.011e+01 6.922e-01 -14.607
private 7.908e+00 7.923e-01 9.982
pub_cap_1000 8.644e-02 1.554e-02 5.562
Correlation of Fixed Effects:
(Intr) year0 yer_sq yer_cb gsp privat
year0 0.122
year_sq 0.021 -0.946
year_cube -0.029 0.884 -0.984
gsp -0.125 -0.039 0.100 -0.124
private -0.506 -0.053 -0.089 0.115 -0.790
pub_cp_1000 0.574 -0.007 0.022 -0.008 -0.184 -0.222
Robustness weights for the residuals:
650 weights are ~= 1. The remaining 166 ones are summarized as
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.240 0.638 0.819 0.768 0.934 0.999
Robustness weights for the random effects:
37 weights are ~= 1. The remaining 11 ones are summarized as
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.451 0.581 0.755 0.725 0.841 0.984
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 (state):
eff: smoothed Huber (k = 1.345, s = 10)
vcp: smoothed Huber, Proposal II (k = 1.345, s = 10) Step 5: Re-Run Model with Outliers Removed (Cook’s Distance > .10)
prod.trimmed <- diagnostics %>%
filter(., .cooksd < .10)
model.trimmed <- lmer(unemp ~ year0 + year_sq + year_cube + gsp + private + pub_cap_1000 + (1|state), REML = FALSE, data = prod.trimmed)Some predictor variables are on very different scales: consider rescalingsummary(model.trimmed)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: unemp ~ year0 + year_sq + year_cube + gsp + private + pub_cap_1000 + (1 | state)
Data: prod.trimmed
AIC BIC logLik deviance df.resid
3034.4 3076.7 -1508.2 3016.4 803
Scaled residuals:
Min 1Q Median 3Q Max
-2.6754 -0.6721 -0.1032 0.5324 3.4192
Random effects:
Groups Name Variance Std.Dev.
state (Intercept) 9.107 3.018
Residual 1.849 1.360
Number of obs: 812, groups: state, 48
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.678e+01 5.157e+00 5.194
year0 -1.779e-02 8.868e-02 -0.201
year_sq 3.859e-02 1.293e-02 2.984
year_cube -1.819e-03 5.316e-04 -3.420
gsp -1.041e+01 7.132e-01 -14.597
private 8.107e+00 8.230e-01 9.850
pub_cap_1000 9.258e-02 1.629e-02 5.685
Correlation of Fixed Effects:
(Intr) year0 yer_sq yer_cb gsp privat
year0 0.134
year_sq 0.019 -0.944
year_cube -0.029 0.882 -0.984
gsp -0.127 -0.045 0.104 -0.127
private -0.517 -0.057 -0.091 0.117 -0.781
pub_cp_1000 0.549 -0.013 0.027 -0.012 -0.165 -0.232
fit warnings:
Some predictor variables are on very different scales: consider rescalingStep 6: Compare Original “Final” Model, Robust Model, and Trimmed Model for Differences (Sensitivity Analysis)
| Model 1 | Model 2 | |
|---|---|---|
| (Intercept) | 27.040 | 26.785 |
| (5.335) | (5.157) | |
| gsp | -10.734 | -10.411 |
| (0.737) | (0.713) | |
| private | 8.409 | 8.107 |
| (0.850) | (0.823) | |
| pub_cap_1000 | 0.091 | 0.093 |
| (0.017) | (0.016) | |
| sd__(Intercept) | 3.115 | 3.018 |
| sd__Observation | 1.412 | 1.360 |
| year_cube | -0.002 | -0.002 |
| (0.001) | (0.001) | |
| year_sq | 0.037 | 0.039 |
| (0.013) | (0.013) | |
| year0 | -0.006 | -0.018 |
| (0.092) | (0.089) | |
| AIC | 3108.9 | 3034.4 |
| BIC | 3151.3 | 3076.7 |
| Log.Lik. | -1545.469 | -1508.212 |
LS0tCnRpdGxlOiAnTXVsdGlsZXZlbCBNb2RlbGluZywgTW9kdWxlIDEwOiBDaGVja2luZyBBc3N1bXB0aW9ucyBmb3IgTUxNcycKYXV0aG9yOiAnRHIuIEJyb2RhJwpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpUaGlzIHdlZWssIHdlIHdpbGwgdXMgdGhlIGBwcm9kdWN0aXZpdHlgIGRhdGFzZXQgZnJvbSBTdGF0YSwgZGV2ZWxvcGVkIGJ5IERyLiBBbGljaWEgTXVubmVsbCBhbmQgcHVibGlzaGVkIGJ5IHRoZSBGZWRlcmFsIFJlc2VydmUgQmFuayBvZiBCb3N0b24gKDE5OTApLgoKIyBMb2FkIGluIE91ciBNVlAgUGFja2FnZXMKYGBge3J9CnN1cHByZXNzUGFja2FnZVN0YXJ0dXBNZXNzYWdlcyhsaWJyYXJ5KHRpZHl2ZXJzZSkpCmxpYnJhcnkobG1lNCkKbGlicmFyeShIbWlzYykKYGBgCgojIFBhcnQgMTogTG9hZGluZyBpbiB0aGUgRGF0YSwgYW5kIHRoZW4gUmVjb2RpbmcgVmFyaWFibGVzCiMjIExvYWQgaW4gdGhlIERhdGEKYGBge3J9Cgpwcm9kdWN0aXZpdHkgPC0gaGF2ZW46OnJlYWRfZHRhKCJwcm9kdWN0aXZpdHkuZHRhIikKCmdsaW1wc2UocHJvZHVjdGl2aXR5KQpgYGAKIyMgVGhlIGBkZXNjcmliZWAgZnVuY3Rpb24gZnJvbSBgSG1pc2NgIGlzIGEgbmljZSBjb2RlYm9vay4uLgpgYGB7cn0KSG1pc2M6OmRlc2NyaWJlKHByb2R1Y3Rpdml0eSkKYGBgCgojIEEgbGl0dGxlIHF1aWNrIGNsZWFuaW5nIC0gQ29udmVydCBSZWdpb24gdG8gRmFjdG9yIFZhcmlhYmxlLCBSZXNjYWxlIFB1YmxpYyBDYXBpdGFsIFN0b2NrLCAiWmVybyBPdXQiIFllYXIgYW5kIENyZWF0ZSBOb25saW5lYXIgVGVybXMKYGBge3J9CnByb2QuY2xlYW4gPC0gcHJvZHVjdGl2aXR5ICU+JQogIG11dGF0ZSguLCAKICAgICAgICAgcmVnaW9uLmZhYyA9IGFzX2ZhY3RvcihyZWdpb24pLAogICAgICAgICB5ZWFyMCA9IHllYXIgLSAxOTcwLAogICAgICAgICB5ZWFyX3NxID0geWVhcjBeMiwKICAgICAgICAgeWVhcl9jdWJlID0geWVhcjBeMywKICAgICAgICAgcHViX2NhcF8xMDAwID0gcHVibGljLzEwMDApCgpnbGltcHNlKHByb2QuY2xlYW4pCmBgYAoKIyBQYXJ0IDI6IEZpbmQgWW91ciAiRmluYWwiIG9yIFByZWZlcnJlZCBNb2RlbAojIyBTdGFydCB3aXRoIGEgbnVsbCBncm93dGggbW9kZWwgc2luY2UgdGhpcyBpcyBsb25naXR1ZGluYWwgZGF0YS4uLgpgYGB7cn0KbW9kZWwubnVsbC5ncm93dGggPC0gbG1lcih1bmVtcCB+IHllYXIwICsgKDF8c3RhdGUpLCBSRU1MID0gRkFMU0UsIGRhdGEgPSBwcm9kLmNsZWFuKQpzdW1tYXJ5KG1vZGVsLm51bGwuZ3Jvd3RoKQpgYGAKCiMjIEFkZCBpbiBRdWFkIEVmZmVjdHMKYGBge3J9Cm1vZGVsLjIgPC0gbG1lcih1bmVtcCB+IHllYXIwICsgeWVhcl9zcSArICgxfHN0YXRlKSwgUkVNTCA9IEZBTFNFLCBkYXRhID0gcHJvZC5jbGVhbikKc3VtbWFyeShtb2RlbC4yKQpgYGAKIyMgQWRkIGluIEN1YmljIEVmZmVjdHMKYGBge3J9Cm1vZGVsLjMgPC0gbG1lcih1bmVtcCB+IHllYXIwICsgeWVhcl9jdWJlICsgKDF8c3RhdGUpLCBSRU1MID0gRkFMU0UsIGRhdGEgPSBwcm9kLmNsZWFuKQpzdW1tYXJ5KG1vZGVsLjMpCmBgYAojIyBCb3RoIFF1YWRyYXRpYyBhbmQgQ3ViaWMgRWZmZWN0cwpgYGB7cn0KbW9kZWwuNCA8LSBsbWVyKHVuZW1wIH4geWVhcjAgKyB5ZWFyX3NxICsgeWVhcl9jdWJlICsgKDF8c3RhdGUpLCBSRU1MID0gRkFMU0UsIGRhdGEgPSBwcm9kLmNsZWFuKQpzdW1tYXJ5KG1vZGVsLjQpCmBgYAoKIyMgTm93IGZvciBzb21lIG1vcmUgcHJlZGljdG9ycy4uLgpgYGB7cn0KbW9kZWwuNSA8LSBsbWVyKHVuZW1wIH4geWVhcjAgKyB5ZWFyX3NxICsgeWVhcl9jdWJlICsgZ3NwICsgcHJpdmF0ZSArIHB1Yl9jYXBfMTAwMCArICgxfHN0YXRlKSwgUkVNTCA9IEZBTFNFLCBkYXRhID0gcHJvZC5jbGVhbikKc3VtbWFyeShtb2RlbC41KQpgYGAKCiMgVXNlIHRoZSBgYXVnbWVudGAgZnVuY3Rpb24gZnJvbSB0aGUgYGJyb29tLm1peGVkYCBQYWNrYWdlIHRvIEdldCBQcmVkaWN0aW9ucywgUmVzaWR1YWxzLCBDb29rJ3MgRGlzdGFuY2VzLCBFdGMuOgpgYGB7cn0KbGlicmFyeShicm9vbS5taXhlZCkKCmRpYWdub3N0aWNzIDwtIGF1Z21lbnQobW9kZWwuNSkKCmBgYAojIFN0ZXAgMTogQXNzZXNzIE5vcm1hbGl0eSBvZiBSZXNpZHVhbHMKIyMgVmlzdWFsbHksIHdpdGggYSBIaXN0b2dyYW0KYGBge3J9CmdncGxvdChkYXRhID0gZGlhZ25vc3RpY3MsIG1hcHBpbmcgPSBhZXMoeCA9IC5yZXNpZCkpICsKICBnZW9tX2hpc3RvZ3JhbShiaW53aWR0aCA9IC4yNSkgKyB0aGVtZV9jbGFzc2ljKCkgKyAKICBsYWJzKHRpdGxlID0gIkhpc3RvZ3JhbSBvZiBSZXNpZHVhbHMgZm9yIFN0YXRlIFVuZW1wbG95bWVudCBNb2RlbCIsCiAgICAgICAgICAgICAgICAgICAgICB4ID0gIlJlc2lkdWFsIFZhbHVlIikgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IGMoLTIuNSwgMi41KSwgbGluZXR5cGUgPSAiZCAgb3R0ZWQiKQpgYGAKIyMgU3RhdGlzdGljYWxseSwgd2l0aCB0aGUgU2hhcGlyby1XaWxrIFRlc3QKYGBge3J9CnNoYXBpcm8udGVzdChkaWFnbm9zdGljcyQucmVzaWQpCmBgYAoKIyBTdGVwIDI6IFVzZSBSZXNpZHVhbHMgdnMuIEZpdHRlZCAoUlZGKSBQbG90IHRvIEFzc2VzcyBIb21vc2tlZGFzdGljaXR5IG9mIEVycm9ycwojIyBCYXNpYyBSVkYgUGxvdCwgV2l0aG91dCBHcm91cHMKYGBge3J9CmdncGxvdChkYXRhID0gZGlhZ25vc3RpY3MsIG1hcHBpbmcgPSBhZXMoeCA9IC5maXR0ZWQsIHkgPSAucmVzaWQpKSArCiAgZ2VvbV9wb2ludCgpICsgbGFicyh0aXRsZSA9ICJSVkYgUGxvdCBmb3IgU3RhdGUgVW5lbXBsb3ltZW50IE1vZGVsIiwKICAgICAgICAgICAgICAgICAgICAgIHggPSAiUHJlZGljdGVkIFZhbHVlLCAlIFVuZW1wbG95bWVudCIsCiAgICAgICAgICAgICAgICAgICAgICB5ID0gIlJlc2lkdWFsIFZhbHVlIikgKyB0aGVtZV9jbGFzc2ljKCkKYGBgCiMjIFNhbWUgVGhpbmcsIEJ1dCBXaXRoIEFkZGVkIENvbG9yIEJ5IFllYXIwCmBgYHtyfQpnZ3Bsb3QoZGF0YSA9IGRpYWdub3N0aWNzLCBtYXBwaW5nID0gYWVzKHggPSAuZml0dGVkLCB5ID0gLnJlc2lkLCBjb2xvciA9IHllYXIwKSkgKwogIGdlb21fcG9pbnQoKSArIGxhYnModGl0bGUgPSAiUlZGIFBsb3QgZm9yIFN0YXRlIFVuZW1wbG95bWVudCBNb2RlbCwgYnkgWWVhciIsCiAgICAgICAgICAgICAgICAgICAgICB4ID0gIlByZWRpY3RlZCBWYWx1ZSwgJSBVbmVtcGxveW1lbnQiLAogICAgICAgICAgICAgICAgICAgICAgeSA9ICJSZXNpZHVhbCBWYWx1ZSIpICsgdGhlbWVfY2xhc3NpYygpCmBgYAojIFN0ZXAgMzogVXNlIENvb2sncyBEaXN0YW5jZSB0byBJZGVudGlmeSBPdXRsaWVycwpgYGB7cn0KZ2dwbG90KGRhdGEgPSBkaWFnbm9zdGljcywgbWFwcGluZyA9IGFlcyh4ID0gLmZpdHRlZCwgeSA9IC5jb29rc2QsIGxhYmVsID0gc3RhdGUpKSArCiAgZ2VvbV9wb2ludCgpICsgZ2VvbV90ZXh0KG51ZGdlX3ggPSAuMjUpICsgdGhlbWVfY2xhc3NpYygpICsgCiAgbGFicyh0aXRsZSA9ICJDb29rJ3MgRGlzdGFuY2UgUGxvdCBmb3IgU3RhdGUgVW5lbXBsb3ltZW50IE1vZGVsIiwKICAgICAgICAgICAgICAgICAgICAgIHggPSAiUHJlZGljdGVkIFZhbHVlLCAlIFVuZW1wbG95bWVudCIsCiAgICAgICAgICAgICAgICAgICAgICB5ID0gIkNvb2sncyBEaXN0YW5jZSIpICsgCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gNC84MTYsIGxpbmV0eXBlID0gImRvdHRlZCIpCiAgCmBgYAoKIyBTdGVwIDQ6IFJlLVJ1biBNb2RlbCB3aXRoIFJvYnVzdCBTdGFuZGFyZCBFcnJvcnMKYGBge3J9CmxpYnJhcnkocm9idXN0bG1tKQptb2RlbC5yb2J1c3QgPC0gcmxtZXIodW5lbXAgfiB5ZWFyMCArIHllYXJfc3EgKyB5ZWFyX2N1YmUgKyBnc3AgKyBwcml2YXRlICsgcHViX2NhcF8xMDAwICsgKDF8c3RhdGUpLCBSRU1MID0gRkFMU0UsIGRhdGEgPSBwcm9kLmNsZWFuKQpzdW1tYXJ5KG1vZGVsLnJvYnVzdCkKYGBgCgojIFN0ZXAgNTogUmUtUnVuIE1vZGVsIHdpdGggT3V0bGllcnMgUmVtb3ZlZCAoQ29vaydzIERpc3RhbmNlID4gLjEwKQpgYGB7cn0KcHJvZC50cmltbWVkIDwtIGRpYWdub3N0aWNzICU+JQogIGZpbHRlciguLCAuY29va3NkIDwgLjEwKQptb2RlbC50cmltbWVkIDwtIGxtZXIodW5lbXAgfiB5ZWFyMCArIHllYXJfc3EgKyB5ZWFyX2N1YmUgKyBnc3AgKyBwcml2YXRlICsgcHViX2NhcF8xMDAwICsgKDF8c3RhdGUpLCBSRU1MID0gRkFMU0UsIGRhdGEgPSBwcm9kLnRyaW1tZWQpCnN1bW1hcnkobW9kZWwudHJpbW1lZCkKCmBgYAojIFN0ZXAgNjogQ29tcGFyZSBPcmlnaW5hbCAiRmluYWwiIE1vZGVsLCBSb2J1c3QgTW9kZWwsIGFuZCBUcmltbWVkIE1vZGVsIGZvciBEaWZmZXJlbmNlcyAoU2Vuc2l0aXZpdHkgQW5hbHlzaXMpCmBgYHtyfQpsaWJyYXJ5KG1vZGVsc3VtbWFyeSkKbW9kZWxzIDwtIGxpc3QobW9kZWwuNSwgbW9kZWwudHJpbW1lZCkKbW9kZWxzdW1tYXJ5KG1vZGVscywgb3V0cHV0ID0gIm1hcmtkb3duIikKYGBgCgo=