HW3

install.packages("multilevel",repos = "http://cran.us.r-project.org")

## 將程式套件安載入 'C:/Users/a0905/AppData/Local/R/win-library/4.2'
## (因為 'lib' 沒有被指定)

## 程式套件 'multilevel' 開啟成功，MD5 和檢查也透過
## 
## 下載的二進位程式套件在
##  C:\Users\a0905\AppData\Local\Temp\RtmpEZL6YQ\downloaded_packages 裡

library(multilevel)

## 載入需要的套件：nlme

## 載入需要的套件：MASS

library(tidyverse)

## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──

## ✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.2      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::collapse() masks nlme::collapse()
## ✖ dplyr::filter()   masks stats::filter()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ dplyr::select()   masks MASS::select()

data(multilevel, "klein2000")

## Warning in data(multilevel, "klein2000"): 沒有 'multilevel' 這個資料集

klein2000 %>% group_by(GRPID) %>% mutate(mneglead=mean(NEGLEAD)) -> dta

summary(m1 <- lme4::lmer(JOBSAT ~ ( PAY | GRPID ), data = dta ))

## Linear mixed model fit by REML ['lmerMod']
## Formula: JOBSAT ~ (PAY | GRPID)
##    Data: dta
## 
## REML criterion at convergence: 3381
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9050 -0.6367 -0.0033  0.6236  3.8665 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  GRPID    (Intercept) 0.8385   0.9157        
##           PAY         1.0682   1.0335   -0.35
##  Residual             4.4165   2.1016        
## Number of obs: 750, groups:  GRPID, 50
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   0.2806     0.1456   1.927

sjPlot::tab_model(m1, show.p=FALSE, show.r2=FALSE)

	JOBSAT
Predictors	Estimates	CI
(Intercept)	0.28	-0.01 – 0.57
Random Effects
σ2	4.42
τ00 GRPID	0.84
τ11 GRPID.PAY	1.07
ρ01 GRPID	-0.35
ICC	0.16
N GRPID	50
Observations	750

summary(m2 <- lme4::lmer(JOBSAT ~ mneglead + (1 | GRPID ), data = dta))

## Linear mixed model fit by REML ['lmerMod']
## Formula: JOBSAT ~ mneglead + (1 | GRPID)
##    Data: dta
## 
## REML criterion at convergence: 3435.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9638 -0.6472 -0.0214  0.6892  3.6020 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  GRPID    (Intercept) 0.3128   0.5593  
##  Residual             5.4741   2.3397  
## Number of obs: 750, groups:  GRPID, 50
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.08414    0.11643   0.723
## mneglead    -1.77696    0.33701  -5.273
## 
## Correlation of Fixed Effects:
##          (Intr)
## mneglead -0.012

sjPlot::tab_model(m2, show.p=FALSE, show.r2=FALSE)

	JOBSAT
Predictors	Estimates	CI
(Intercept)	0.08	-0.14 – 0.31
mneglead	-1.78	-2.44 – -1.12
Random Effects
σ2	5.47
τ00 GRPID	0.31
ICC	0.05
N GRPID	50
Observations	750

anova(m1, m2)

## refitting model(s) with ML (instead of REML)

## Data: dta
## Models:
## m2: JOBSAT ~ mneglead + (1 | GRPID)
## m1: JOBSAT ~ (PAY | GRPID)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m2    4 3440.3 3458.8 -1716.2   3432.3                         
## m1    5 3389.0 3412.1 -1689.5   3379.0 53.386  1  2.741e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1