Please perform multilevel analyses using below models in hierarchical forms. Use maximum likelihood (ML) estimation and unconstructured covariance type. Please write down the combined form and all estimated values for γ (fixed effect), σ2, τ (random effects). You don’t have to interpret them.
# Import the dataset
library(haven)
HM6 <- read_sav("hsbdataset_HM1.sav")
View(HM6)
# Preview the structure of the data
str(HM6)## Classes 'tbl_df', 'tbl' and 'data.frame': 7185 obs. of 14 variables:
## $ school : num 1224 1224 1224 1224 1224 ...
## ..- attr(*, "format.spss")= chr "F10.0"
## $ student : num 1 2 3 4 5 6 7 8 9 10 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ cons : num 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ minority: num 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ female : num 1 1 0 0 0 0 1 0 1 0 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ ses : num -1.528 -0.588 -0.528 -0.668 -0.158 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ meanses : num -0.428 -0.428 -0.428 -0.428 -0.428 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ cses : num -1.1 -0.16 -0.1 -0.24 0.27 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ mathach : num 5.88 19.71 20.35 8.78 17.9 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ size : num 842 842 842 842 842 842 842 842 842 842 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ sector : num 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ pracad : num 0.35 0.35 0.35 0.35 0.35 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ disclim : num 1.6 1.6 1.6 1.6 1.6 ...
## ..- attr(*, "format.spss")= chr "F9.2"
## $ himinty : num 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "format.spss")= chr "F9.2"# Preview the descriptives of the data
summary(HM6)## school student cons minority
## Min. :1224 Min. : 1.00 Min. :1 Min. :0.0000
## 1st Qu.:3020 1st Qu.:12.00 1st Qu.:1 1st Qu.:0.0000
## Median :5192 Median :23.00 Median :1 Median :0.0000
## Mean :5278 Mean :24.51 Mean :1 Mean :0.2747
## 3rd Qu.:7342 3rd Qu.:36.00 3rd Qu.:1 3rd Qu.:1.0000
## Max. :9586 Max. :67.00 Max. :1 Max. :1.0000
## female ses meanses
## Min. :0.0000 Min. :-3.758000 Min. :-1.188000
## 1st Qu.:0.0000 1st Qu.:-0.538000 1st Qu.:-0.317000
## Median :1.0000 Median : 0.002000 Median : 0.038000
## Mean :0.5282 Mean : 0.000143 Mean : 0.006138
## 3rd Qu.:1.0000 3rd Qu.: 0.602000 3rd Qu.: 0.333000
## Max. :1.0000 Max. : 2.692000 Max. : 0.831000
## cses mathach size sector
## Min. :-3.657000 Min. :-2.832 Min. : 100 Min. :0.0000
## 1st Qu.:-0.454000 1st Qu.: 7.275 1st Qu.: 565 1st Qu.:0.0000
## Median : 0.010000 Median :13.131 Median :1016 Median :0.0000
## Mean :-0.005995 Mean :12.748 Mean :1057 Mean :0.4931
## 3rd Qu.: 0.463000 3rd Qu.:18.317 3rd Qu.:1436 3rd Qu.:1.0000
## Max. : 2.850000 Max. :24.993 Max. :2713 Max. :1.0000
## pracad disclim himinty
## Min. :0.0000 Min. :-2.4160 Min. :0.00
## 1st Qu.:0.3200 1st Qu.:-0.8170 1st Qu.:0.00
## Median :0.5300 Median :-0.2310 Median :0.00
## Mean :0.5345 Mean :-0.1319 Mean :0.28
## 3rd Qu.:0.7000 3rd Qu.: 0.4600 3rd Qu.:1.00
## Max. :1.0000 Max. : 2.7560 Max. :1.00# Prepare the ICC function for future use
ICC <- function(out) {
varests <- as.numeric(VarCorr(out)[1:nrow(VarCorr(out))])
return(paste("ICC=",varests[1]/ (varests[1]+varests[-1])))
}# For model 1, we need to build a random intercept and slope model with one interaction terms.
library(nlme)
Model_1 <- lme(mathach~meanses + sector + cses + cses:himinty + female, random=~cses|school,data=HM6,method="ML")
summary(Model_1)## Linear mixed-effects model fit by maximum likelihood
## Data: HM6
## AIC BIC logLik
## 46499.86 46568.65 -23239.93
##
## Random effects:
## Formula: ~cses | school
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 1.4508784 (Intr)
## cses 0.7423997 0.123
## Residual 6.0420327
##
## Fixed effects: mathach ~ meanses + sector + cses + cses:himinty + female
## Value Std.Error DF t-value p-value
## (Intercept) 12.681442 0.2074586 7022 61.12758 0.0000
## meanses 5.184546 0.3528852 157 14.69188 0.0000
## sector 1.329812 0.2923686 157 4.54841 0.0000
## cses 2.337350 0.1484546 7022 15.74454 0.0000
## female -1.182440 0.1617618 7022 -7.30976 0.0000
## cses:himinty -0.614294 0.2712710 7022 -2.26450 0.0236
## Correlation:
## (Intr) meanss sector cses female
## meanses 0.207
## sector -0.635 -0.356
## cses 0.014 0.012 -0.005
## female -0.400 0.044 -0.011 0.046
## cses:himinty 0.000 -0.024 0.011 -0.546 -0.017
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.18176524 -0.73024820 0.02390801 0.75313970 2.90449106
##
## Number of Observations: 7185
## Number of Groups: 160anova(Model_1)## numDF denDF F-value p-value
## (Intercept) 1 7022 8511.586 <.0001
## meanses 1 157 315.871 <.0001
## sector 1 157 20.016 <.0001
## cses 1 7022 311.668 <.0001
## female 1 7022 54.000 <.0001
## cses:himinty 1 7022 5.128 0.0236logLik(Model_1)*-2 ## 'log Lik.' 46479.86 (df=10)VarCorr(Model_1)## school = pdLogChol(cses)
## Variance StdDev Corr
## (Intercept) 2.1050482 1.4508784 (Intr)
## cses 0.5511573 0.7423997 0.123
## Residual 36.5061592 6.0420327Model in combined form:
Y_ij=(γ_00+ γ_01 M〖EAN SEAS〗_J+ γ_02 〖SECTOR〗_j+ γ_10 〖female〗_ij+ γ_11 female〖MEAN SES〗_j+ γ_12 female〖SECTOR〗_j)+ (u_0j+〖ufemale〗_1j+r_ij)
Estimates that could be found in the previous outcomes:
# For model 1, we need to build a random intercept and slope model with two interaction terms.
Model_2 <- lme(mathach~meanses + sector + female + meanses*female + sector:female , random=~female|school,data=HM6,method="ML")
summary(Model_2)## Linear mixed-effects model fit by maximum likelihood
## Data: HM6
## AIC BIC logLik
## 46890.61 46959.41 -23435.31
##
## Random effects:
## Formula: ~female | school
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 1.6705983 (Intr)
## female 0.8005715 -0.699
## Residual 6.2242449
##
## Fixed effects: mathach ~ meanses + sector + female + meanses * female + sector:female
## Value Std.Error DF t-value p-value
## (Intercept) 12.853537 0.2416952 7022 53.18077 0.0000
## meanses 5.220786 0.4841625 157 10.78313 0.0000
## sector 1.082889 0.3979865 157 2.72092 0.0072
## female -1.472157 0.2360288 7022 -6.23719 0.0000
## meanses:female -0.053704 0.5011990 7022 -0.10715 0.9147
## sector:female 0.274808 0.4139581 7022 0.66385 0.5068
## Correlation:
## (Intr) meanss sector female mnss:f
## meanses 0.255
## sector -0.675 -0.419
## female -0.639 -0.193 0.449
## meanses:female -0.178 -0.702 0.326 0.291
## sector:female 0.420 0.330 -0.695 -0.652 -0.449
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.23125355 -0.74522523 0.03426301 0.76539543 2.69003656
##
## Number of Observations: 7185
## Number of Groups: 160anova(Model_2)## numDF denDF F-value p-value
## (Intercept) 1 7022 8707.631 <.0001
## meanses 1 157 325.063 <.0001
## sector 1 157 19.024 <.0001
## female 1 7022 58.593 <.0001
## meanses:female 1 7022 0.046 0.8311
## sector:female 1 7022 0.441 0.5068logLik(Model_2)*-2 ## 'log Lik.' 46870.61 (df=10)VarCorr(Model_2)## school = pdLogChol(female)
## Variance StdDev Corr
## (Intercept) 2.7908985 1.6705983 (Intr)
## female 0.6409147 0.8005715 -0.699
## Residual 38.7412244 6.2242449Model in combined form:
Y_ij=(γ_00+ γ_01 M〖EAN SEAS〗_J + γ_02〖SECTOR〗_j + γ_10〖CSES〗_i j+ γ_11 himinty〖CES〗_ij+ γ_20〖Female〗_ij )+ (u_0j +〖uCSES〗_1j + r_ij)
Estimates that could be found in the previous outcomes: