Thực hành mô hình đa tầng
DỰng mô hình đa tầng với ý tưởng mỗi con chuột có bắt đầu khác nhau (alpha0=A0+a) trong đó A0 là hằng số a là sai số VÀ thay đổi theo thời gian là TIME = Time(hằng số và Time sai số) công Thức tổng để thay 2 cái trên vào là Glucose = Hằng số khởi đầu + Biến thiên thời gian. Ta có:
library(lme4)
## Loading required package: Matrix
attach(Glucose)
m = glmer(glucose~time + (1 + time|id), data=Glucose)
## Warning in glmer(glucose ~ time + (1 + time | id), data = Glucose): calling
## glmer() with family=gaussian (identity link) as a shortcut to lmer() is
## deprecated; please call lmer() directly
## boundary (singular) fit: see ?isSingular
summary(m)
## Linear mixed model fit by REML ['lmerMod']
## Formula: glucose ~ time + (1 + time | id)
## Data: Glucose
##
## REML criterion at convergence: 165.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6169 -0.4435 -0.0086 0.4548 3.2935
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.6900322 0.83068
## time 0.0009616 0.03101 1.00
## Residual 0.2536389 0.50363
## Number of obs: 76, groups: id, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.0481 0.2177 27.79
## time -0.4699 0.0397 -11.84
##
## Correlation of Fixed Effects:
## (Intr)
## time -0.240
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Công thức tổng như cũ, 2 công thức tầng 2 chứa biến treatment trở thành (ví dụ: TIME= hằng số + Time.treat + sai số; tương tự với khởi đầu)
m2 = glmer( glucose~treatment+time+treatment*time + (1+time|id), data=Glucose)
## Warning in glmer(glucose ~ treatment + time + treatment * time + (1 + time | :
## calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is
## deprecated; please call lmer() directly
## boundary (singular) fit: see ?isSingular
summary(m2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: glucose ~ treatment + time + treatment * time + (1 + time | id)
## Data: Glucose
##
## REML criterion at convergence: 166.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.60331 -0.41198 0.00053 0.48474 2.99267
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.7680707 0.87640
## time 0.0001578 0.01256 1.00
## Residual 0.2534072 0.50340
## Number of obs: 76, groups: id, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.96539 0.29442 20.262
## treatment 0.18082 0.39045 0.463
## time -0.41860 0.05327 -7.859
## treatment:time -0.11451 0.07910 -1.448
##
## Correlation of Fixed Effects:
## (Intr) trtmnt time
## treatment -0.637
## time -0.403 0.357
## treatmnt:tm 0.269 -0.437 -0.675
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Lúc này có thể dùng ggplot vẽ hình nhẹ về việc fill bạn treatment rồi group by bạn id để thấy được hiểu quả nghiên cứu, có thể box, có thể spaghetti
Mô hình trên do ảnh hưởng của sai số theo thời gian quá nhỏ nên hủy luôn biến số time trong sai số, làm lại thành:
m3 = glmer( glucose~treatment+time+treatment:time + (1|id), data=Glucose)
## Warning in glmer(glucose ~ treatment + time + treatment:time + (1 | id), :
## calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is
## deprecated; please call lmer() directly
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: glucose ~ treatment + time + treatment:time + (1 | id)
## Data: Glucose
##
## REML criterion at convergence: 166.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.60263 -0.43223 -0.00001 0.53418 2.95623
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.8148 0.9027
## Residual 0.2542 0.5042
## Number of obs: 76, groups: id, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.96580 0.30178 19.769
## treatment 0.18047 0.40102 0.450
## time -0.41616 0.05342 -7.790
## treatment:time -0.12016 0.07925 -1.516
##
## Correlation of Fixed Effects:
## (Intr) trtmnt time
## treatment -0.638
## time -0.453 0.384
## treatmnt:tm 0.295 -0.476 -0.678
Một data mới. Nghiên cứu về yếu tố ảnh hưởng điểm học Nhóm nghiên cứu muốn đánh giá xem mức độ thường xuyên làm bài tập ở nhà ảnh hưởng như thế nào đến điểm thi toán. Thông tin thu nhận từ 260 học sinh của 10 trường. Giới thiệu tập tin số liệu imm10.csv: - id: mã số học sinh tham gia nghiên cứu - schnum: mã số trường - math: điểm thi toán - homework: mức độ thường xuyên làm bài tập ở nhà, tính bằng điểm dao động từ 0 đến 5 - white: người da trắng - sex: giới tính (1= Nam, 2= Nữ) - Các yếu tố ảnh hưởng khác: ses, sex, sctype, sctr, scsize, urban, region (không sử dụng trong bài tập này)
imm10 = read.csv("C:/Users/Paul Vo/Desktop/Textbook/TDTU Datasets for 2020 Workshop/imm10.csv")
head(imm10)
## school id ses meanses homework white parented public ratio percmin math
## 1 7472 3 -0.13 -0.4826087 1 1 2 1 19 0 48
## 2 7472 8 -0.39 -0.4826087 0 1 2 1 19 0 48
## 3 7472 13 -0.80 -0.4826087 0 1 2 1 19 0 53
## 4 7472 17 -0.72 -0.4826087 1 1 2 1 19 0 42
## 5 7472 27 -0.74 -0.4826087 2 1 2 1 19 0 43
## 6 7472 28 -0.58 -0.4826087 1 1 2 1 19 0 57
## sex race sctype cstr scsize urban region schnum
## 1 2 4 1 2 3 2 2 1
## 2 1 4 1 2 3 2 2 1
## 3 1 4 1 2 3 2 2 1
## 4 1 4 1 2 3 2 2 1
## 5 2 4 1 2 3 2 2 1
## 6 2 4 1 2 3 2 2 1
Mô hình khi Bạn muốn đánh giá ảnh hưởng của làm bài tập ở nhà lên điểm thi toán dưới giả định là học sinhcó làm bài tập ảnh hưởng thế nào nếu xét cùng với trường
library(lme4)
mi = lmer(math~homework+(1+homework|school), data = imm10)
summary(mi)
## Linear mixed model fit by REML ['lmerMod']
## Formula: math ~ homework + (1 + homework | school)
## Data: imm10
##
## REML criterion at convergence: 1764
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5111 -0.5357 0.0175 0.6121 2.5708
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## school (Intercept) 69.32 8.326
## homework 22.46 4.739 -0.81
## Residual 43.07 6.563
## Number of obs: 260, groups: school, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 44.771 2.744 16.316
## homework 2.040 1.554 1.313
##
## Correlation of Fixed Effects:
## (Intr)
## homework -0.804
Bạn muốn đánh giá ảnh hưởng của làm bài tập ở nhà lên odds học sinh giỏi (định nghĩa làm điểm toán> 65) dưới giả định là học sinh tham gia nghiên cứu độc lập với nhau (~ không có khác biệt giữa học sinh các trường): 4.1 Làm bài tập ở nhà ảnh hưởng như thế nào đến odds học sinh giỏi?
imm10$good= ifelse(imm10$math>65, 1, 0)
mi2= glm(good~ homework, family= binomial, data= imm10)
summary(mi2)
##
## Call:
## glm(formula = good ~ homework, family = binomial, data = imm10)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8182 -0.3986 -0.2719 -0.2719 2.5755
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.0662 0.4901 -8.297 < 2e-16 ***
## homework 0.7867 0.1334 5.898 3.68e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 197.86 on 259 degrees of freedom
## Residual deviance: 154.41 on 258 degrees of freedom
## AIC: 158.41
##
## Number of Fisher Scoring iterations: 5
Lặp lại câu 4 dưới giả định là tình hình học toán và làm bài tập ở nhà khác nhau giữa các trường.
mi3= lmer(good~ homework+ (1+homework | school), data= imm10)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00657263 (tol = 0.002, component 1)
summary(mi3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: good ~ homework + (1 + homework | school)
## Data: imm10
##
## REML criterion at convergence: 107.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5377 -0.2971 -0.0914 -0.0118 3.2444
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## school (Intercept) 0.0001061 0.01030
## homework 0.0023573 0.04855 -1.00
## Residual 0.0812302 0.28501
## Number of obs: 260, groups: school, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.02138 0.03152 -0.678
## homework 0.05475 0.02219 2.467
##
## Correlation of Fixed Effects:
## (Intr)
## homework -0.635
## convergence code: 0
## Model failed to converge with max|grad| = 0.00657263 (tol = 0.002, component 1)