BÀI THỰC HÀNH SỐ 4 - MÔ HÌNH ĐA TẦNG
Phương Anh
8/1/2020
Q1. Nghiên cứu về hiệu quả thuốc lên nồng độ đường huyết của chuột
– Nhóm nghiên cứu muốn đánh giá xem liệu sử dụng thuốc A có hiệu quả làm giảm đường huyết không. Nghiên cứu thực hiện trên 19 chuột, trong đó 10 chuột không dùng thuốc và 9 chuột được dùng thuốc. Nồng độ đường huyết của chuột được đo ở thời điểm trước và 3 thời điểm sau dùng thuốc.
– Giới thiệu tập tin số liệu glucose.csv:
* id: mã số của chuột
* time: thời điểm đo nồng độ đường huyết (trước khi dùng thuốc (T0), 2 (T2), 3 (T3) và 4 giờ (T4) sau dùng thuốc).
* treatment: can thiệp (1= có dùng thuốc, 0= không dùng thuốc)
#Đọc dữ liệu
t = "F:\\NCKH - PHAN TICH THONG KE\\BAI GIANG CAC KHOA PT DL\\Kho hoc - Du lieu da tang SG1.2020\\Data for practice\\glucose.csv"
gl = read.csv(t, header = T)
head(gl)## id treatment time glucose
## 1 1 1 0 5.9
## 2 1 1 2 3.9
## 3 1 1 3 3.9
## 4 1 1 4 3.6
## 5 2 1 0 5.3
## 6 2 1 2 4.71.1 Nồng độ đường huyết của các chuột tham gia nghiên cứu thay đổi qua từng thời điểm như thế nào?
library(ggplot2)
#Vẽ graph về glucose theo cá thể cho tất cả (gồm nhóm có và không can thiệp)
gl$treatment = as.factor(gl$treatment)
ggplot(data = gl, aes(x=time, y = glucose, group = id, col = treatment)) +
geom_point(size = 3) +
geom_line(col = "black")
#Vẽ graph về glucose theo cá thể cho 2 nhóm riêng (có và không can thiệp)
library(gridExtra)
ggplot(data= gl) +
geom_point(aes(x= time,y= glucose),
color="blue", size=3)+
geom_line(aes(x= time, y= glucose, group= id))+
facet_grid(~ treatment) +
stat_smooth(aes(x= time, y= glucose))## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.02## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 3.02## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 1.6401e-016## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0804## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.02## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 3.02## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 1.6401e-016## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0804## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.02## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 3.02## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 2.2134e-016## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0804## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.02## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 3.02## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 2.2134e-016## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0804
1.2 Đường huyết khi chưa uống thuốc của chuột ở 2 nhóm có khác nhau không? Làm sao để xác định sự khác nhau này không phải do ngẫu nhiên?
#Tách dữ liệu
##Cách 1: phương pháp đơn giản/thủ công
head(gl)## id treatment time glucose
## 1 1 1 0 5.9
## 2 1 1 2 3.9
## 3 1 1 3 3.9
## 4 1 1 4 3.6
## 5 2 1 0 5.3
## 6 2 1 2 4.7t0 = subset(gl, time == 0)
t2 = subset(gl, time == 2)
t0$glucose0 = t0$glucose
t2$glucose2 = t2$glucose
t0 = t0[, c("id", "treatment", "glucose0")]
t2 = t2[, c("id", "treatment", "glucose2")]
t = merge(t0,t2, by = c("id", "treatment"))
head(t)## id treatment glucose0 glucose2
## 1 1 1 5.9 3.9
## 2 10 0 6.2 5.3
## 3 11 0 6.9 5.6
## 4 12 0 5.6 4.7
## 5 13 0 5.1 3.9
## 6 14 0 5.7 4.7head(t)## id treatment glucose0 glucose2
## 1 1 1 5.9 3.9
## 2 10 0 6.2 5.3
## 3 11 0 6.9 5.6
## 4 12 0 5.6 4.7
## 5 13 0 5.1 3.9
## 6 14 0 5.7 4.7head(t0)## id treatment glucose0
## 1 1 1 5.9
## 5 2 1 5.3
## 9 3 1 4.6
## 13 4 1 6.0
## 17 5 1 6.0
## 21 6 1 6.4head(t2)## id treatment glucose2
## 2 1 1 3.9
## 6 2 1 4.7
## 10 3 1 3.7
## 14 4 1 5.4
## 18 5 1 5.4
## 22 6 1 4.7##Cách 2 : PP hiện đại đang được sử dụng, dạng level cao
###Lấy 1 phần dữ liệu gl có time = 0 hoặc time = 2
library(tidyverse)## -- Attaching packages --------------------------------------------------------------------- tidyverse 1.3.0 --## v tibble 2.1.3 v dplyr 0.8.3
## v tidyr 1.0.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.4.0
## v purrr 0.3.3## -- Conflicts ------------------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::combine() masks gridExtra::combine()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()g0 = filter(.data = gl, time ==0 )
g2 = filter(.data = gl, time == 2)
g02 = filter(.data = gl, time ==0 | time == 2)
###Chuyển dữ liệu từ dạng dọc sang dạng ngang với biến time
library(reshape2)##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smithsg3 = dcast(g02, id + treatment ~ time)## Using glucose as value column: use value.var to override.head(g3)## id treatment 0 2
## 1 1 1 5.9 3.9
## 2 2 1 5.3 4.7
## 3 3 1 4.6 3.7
## 4 4 1 6.0 5.4
## 5 5 1 6.0 5.4
## 6 6 1 6.4 4.7#vẽ boxplot giữa theo nhóm treatment tại thời điểm t0
p = ggplot(data = t0, aes(x=factor(treatment), y = glucose0, fill = treatment))
p = p + geom_boxplot()
p
p = ggplot(data = g0, aes(x=factor(treatment), y = glucose, fill = treatment))
p = p + geom_boxplot()
p
#vẽ boxplot giữa theo nhóm treatment tại thời điểm t0, thêm tất cả dữ liệu
p = ggplot(data = t0, aes(x=factor(treatment), y = glucose0, fill = treatment))
p = p + geom_boxplot() + geom_jitter(alpha = 0.2, size = 2)
p
p = ggplot(data = g0, aes(x=factor(treatment), y = glucose, fill = treatment))
p = p + geom_boxplot() + geom_jitter(alpha = 0.2, size = 2)
p
#vẽ boxplot giữa theo nhóm treatment tại thời điểm t0 và t2
p = ggplot(data = g2, aes(x=factor(treatment), y = glucose, fill = treatment))
p = p + geom_boxplot() + geom_jitter(alpha = 0.2, size = 2)
p
#Kiểm tra khác biệt giữa 2 nhóm điều trị và không điều trị tại thời điểm ban đầu t0
library(coin)## Loading required package: survivaloneway_test(g0$glucose ~ factor(g0$treatment))##
## Asymptotic Two-Sample Fisher-Pitman Permutation Test
##
## data: g0$glucose by factor(g0$treatment) (0, 1)
## Z = 0.36664, p-value = 0.7139
## alternative hypothesis: true mu is not equal to 0median_test(g0$glucose ~ factor(g0$treatment), conf.int = TRUE)##
## Asymptotic Two-Sample Brown-Mood Median Test
##
## data: g0$glucose by factor(g0$treatment) (0, 1)
## Z = 0.7151, p-value = 0.4745
## alternative hypothesis: true mu is not equal to 0
## 95 percent confidence interval:
## -0.8 1.8
## sample estimates:
## difference in location
## 0.21.3 Bạn dự định làm thế nào để đánh giá hiệu quả của thuốc lên nồng độ đường máu?
1.3.1 Cách 1: giả thiết là bạn sử dụng mức độ khác biệt giữa 2 nhóm về nồng độ đường máu ở thời điểm 2 giờ sau dùng thuốc. Thực hiện trên R và nhận xét kết quả.
#Lọc dữ liệu với time== 0 và time== 2
g02 = filter(.data = gl, time ==0 | time == 2)
#Biến đổi dữ liệu từ dòng sang cột bằng package reshape2 và hàm dcast
g3 = dcast(g02, id + treatment ~ time)## Using glucose as value column: use value.var to override.head(g3)## id treatment 0 2
## 1 1 1 5.9 3.9
## 2 2 1 5.3 4.7
## 3 3 1 4.6 3.7
## 4 4 1 6.0 5.4
## 5 5 1 6.0 5.4
## 6 6 1 6.4 4.7#Đổi tên biến số 0 thành glucose0 và 2 thành glucose2
names(g3)[names(g3) == "0"] = "glucose0"
names(g3)[names(g3) == "2"] = "glucose2"
head(g3)## id treatment glucose0 glucose2
## 1 1 1 5.9 3.9
## 2 2 1 5.3 4.7
## 3 3 1 4.6 3.7
## 4 4 1 6.0 5.4
## 5 5 1 6.0 5.4
## 6 6 1 6.4 4.7#Kiểm định sự khác biệt về mức đường huyết ở thời điểm 2 giờ ở nhóm điều trị và không điều trị
oneway_test(g3$glucose2 ~ factor(g0$treatment))##
## Asymptotic Two-Sample Fisher-Pitman Permutation Test
##
## data: g3$glucose2 by factor(g0$treatment) (0, 1)
## Z = 0.93739, p-value = 0.3486
## alternative hypothesis: true mu is not equal to 0median_test(g3$glucose2 ~ factor(g0$treatment), conf.int = TRUE)##
## Asymptotic Two-Sample Brown-Mood Median Test
##
## data: g3$glucose2 by factor(g0$treatment) (0, 1)
## Z = 0.7151, p-value = 0.4745
## alternative hypothesis: true mu is not equal to 0
## 95 percent confidence interval:
## -0.7 1.7
## sample estimates:
## difference in location
## 0.31.3.2 Cách 2: giá thiết là bạn đánh giá mức độ khác biệt giữa 2 nhóm về thay đổi đường huyết cho đến thời điểm 2 giờ sau dùng thuốc (thay đổi đường huyết= nồng độ đường huyết thời điểm 2 giờ- nồng độ đường huyết khi chưa dùng thuốc). Thực hiện trên R và nhận xét kết quả.
#Đặt biến thay đổi đường huyết là diff= glucose2 - glucose0
g3$diff = g3$glucose2 - g3$glucose0
head(g3)## id treatment glucose0 glucose2 diff
## 1 1 1 5.9 3.9 -2.0
## 2 2 1 5.3 4.7 -0.6
## 3 3 1 4.6 3.7 -0.9
## 4 4 1 6.0 5.4 -0.6
## 5 5 1 6.0 5.4 -0.6
## 6 6 1 6.4 4.7 -1.7#Kiểm định sự khác biệt về thay đổi đường huyết ở thời điểm 2 giờ ở nhóm điều trị và không điều trị
oneway_test(g3$diff ~ factor(g3$treatment))##
## Asymptotic Two-Sample Fisher-Pitman Permutation Test
##
## data: g3$diff by factor(g3$treatment) (0, 1)
## Z = 0.51813, p-value = 0.6044
## alternative hypothesis: true mu is not equal to 01.4 Bạn sử dụng mô hình phân tích ảnh hưởng hỗn hợp để trả lời câu hỏi nghiên cứu
1.4.1 Bạn phát biểu câu hỏi nghiên cứu chính xác bạn dự định trả lời bằng mô hình ảnh hưởng hỗn hợp. Câu hỏi này khác như thế nào với phần 1.3 trên.
1.4.2 Bạn hãy xây dựng mô hình đánh giá hiệu quả của thuốc trên nồng độ đường huyết của chuột i bằng những thông tin sau:
- yij: nồng độ đường huyết của chuột i ở thời điểm j
- A0: nồng độ trung bình của đường huyết trước uống thuốc của chuột trong nhóm chứng
- A1: khác biệt về nồng độ đường huyết trước uống thuốc giữa 2 nhóm
- B0: tốc độ thay đổi trung bình của đường huyết của chuột trong nhóm chứng
- B1: khác biệt về tốc độ thay đổi nồng độ đường huyết giữa 2 nhóm
- u: dao động nồng độ đường huyết trước uống thuốc giữa các chuột
- v: dao động tốc độ thay đổi đường huyết giữa các chuột
- eij: giá trị nhiễu (≠ nồng độ đường huyết giữa giá trị đo được và giá trị dự báo từ mô hình) của chuột i.
- treatment (điểu trị); time (thời điểm đo đường huyết j)
Mô hình Unconditional Model - xét ảnh hưởng của time, không xét ảnh hưởng của treatment :
* Yij = Alphai + Betai + eij * Alphai = A0 + random u * Betai = B0 + random v * Yij = (A0 + u) + (B0 + v).timeij + eij * Yij = (A0 + B0.timeij) + (u + v.timeij + eij)
Mô hình Conditional Model : xét ảnh hưởng của time, xét thêm ảnh hưởng của treatment
* Alphai = A0 + A1.Treatment + ui
* Betai = B0 + B1.Treatment + vi
* Yij = (A0 + A1.Treat + u) + (B0 + B1.Treat + v).Timeij + eij
* Yij = (A0 + A1.Treat + B0.Time + B1.Treat.Time) + (u + v.Timeij + eij)
1.4.3 Thực hiện mô hình phân tích ảnh hưởng hỗn hợp bằng R
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpackm = lmer(glucose ~ time + (1 + time|id), data = gl)## boundary (singular) fit: see ?isSingularsummary(m)## Linear mixed model fit by REML ['lmerMod']
## Formula: glucose ~ time + (1 + time | id)
## Data: gl
##
## 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 ?isSingularm2 = lmer(glucose ~ treatment + time + treatment*time + (1 + time|id), data = gl)## boundary (singular) fit: see ?isSingularsummary(m2)## Linear mixed model fit by REML ['lmerMod']
## Formula: glucose ~ treatment + time + treatment * time + (1 + time | id)
## Data: gl
##
## 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
## treatment1 0.18082 0.39045 0.463
## time -0.41860 0.05327 -7.859
## treatment1:time -0.11451 0.07910 -1.448
##
## Correlation of Fixed Effects:
## (Intr) trtmn1 time
## treatment1 -0.637
## time -0.403 0.357
## tretmnt1:tm 0.269 -0.437 -0.675
## convergence code: 0
## boundary (singular) fit: see ?isSingular1.4.4 Viết lại mô hình 1.4.2 bằng các ước số từ kết quả mô hình
Glucose = 5.97 + 0.18.Treatment - 0.42.Time - 0.11.Treatment.Time
- Nhóm Control : 5.97 - 0.42.Time
- Nhóm Treatment: 6.15 - 0.56.Time
+ A0 = 5.96 : nồng độ trung bình của đường huyết trước uống thuốc của chuột trong nhóm chứng, là nồng độ glucose lúc ban đầu của nhóm chứng
+ A1 = 0.18 : khác biệt về nồng độ TB của glucose (đường huyết) ban đầu (trước uống thuốc) giữa 2 nhóm
+ B0 = -0.42 : tốc độ thay đổi trung bình của đường huyết của chuột trong nhóm chứng
+ B1 = -0.11 : khác biệt về tốc độ thay đổi nồng độ đường huyết giữa 2 nhóm
+ u = 0.77 : dao động nồng độ đường huyết trước uống thuốc giữa các chuột, phản ánh khác biệt giữa các giá trị ban đầu
+ v = 0.00016 : dao động tốc độ thay đổi đường huyết giữa các chuột, phản ánh khác biệt giữa các slope của các con chuột –> Rất nhỏ, khác biệt không đáng kể. Điều này chứng tỏ mặc dù có sự dao động giữa 2 nhóm nhưng mức độ khác biệt giữa các con chuột không đáng kể –> có thể mô hình lại bỏ v
+ eij = 0.25 : giá trị nhiễu (≠ nồng độ đường huyết giữa giá trị đo được và giá trị dự báo từ mô hình) của chuột i.
# Mô hình lại với việc bỏ qua khác biệt slope giữa các con chuột (v)
m2v = lmer(glucose ~ treatment + time + treatment*time + (1|id), data = gl)
summary(m2v)## Linear mixed model fit by REML ['lmerMod']
## Formula: glucose ~ treatment + time + treatment * time + (1 | id)
## Data: gl
##
## 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
## treatment1 0.18047 0.40102 0.450
## time -0.41616 0.05342 -7.790
## treatment1:time -0.12016 0.07925 -1.516
##
## Correlation of Fixed Effects:
## (Intr) trtmn1 time
## treatment1 -0.638
## time -0.453 0.384
## tretmnt1:tm 0.295 -0.476 -0.678Q2. Nghiên cứu Orthodont
– Nhằm xác định xem có sự khác biệt về tốc độ phát triển giữa nam và nữ, nhóm nghiên cứu theo dõi phát triển của 17 trẻ em (16 nam và 11 nữ) từ 8 đến 14 tuổi. Mỗi 2 năm nghiên cứu viên đo khoảng cách giữa tuyến yên (pituitary) và khe chân bướm hàm (pterygomaxillary) trên film x-quang (distance). Tập tin Orthodont (có sẵn trong gói nlme) bao gồm 108 dòng và 4 cột. Bốn biến số gồm:
+ distance: khoảng cách (mm) giữa tuyến yên và khe chân bướm hàm đo được trên film x-quang.
+ age: tuổi (năm)
+ Subject: mã số của người tham gia (M01-M16: nam; F01-F13: nữ).
+ Sex: giới tính (Male: nam; Female: nữ)
2.1 Nhận xét sơ lược về khác biệt giữa nam và nữ về:
(i) khoảng cách distance ở thời điểm bắt đầu tham gia nghiên cứu
(ii) thay đổi distance theo thời gian
(iii) các dạng thay đổi bất thường.
#ĐỌc dữ liệu
data("Orthodont", package="nlme")
head(Orthodont)## Grouped Data: distance ~ age | Subject
## distance age Subject Sex
## 1 26.0 8 M01 Male
## 2 25.0 10 M01 Male
## 3 29.0 12 M01 Male
## 4 31.0 14 M01 Male
## 5 21.5 8 M02 Male
## 6 22.5 10 M02 Malestr(Orthodont)## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 108 obs. of 4 variables:
## $ distance: num 26 25 29 31 21.5 22.5 23 26.5 23 22.5 ...
## $ age : num 8 10 12 14 8 10 12 14 8 10 ...
## $ Subject : Ord.factor w/ 27 levels "M16"<"M05"<"M02"<..: 15 15 15 15 3 3 3 3 7 7 ...
## $ Sex : Factor w/ 2 levels "Male","Female": 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "outer")=Class 'formula' language ~Sex
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## - attr(*, "formula")=Class 'formula' language distance ~ age | Subject
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## - attr(*, "labels")=List of 2
## ..$ x: chr "Age"
## ..$ y: chr "Distance from pituitary to pterygomaxillary fissure"
## - attr(*, "units")=List of 2
## ..$ x: chr "(yr)"
## ..$ y: chr "(mm)"
## - attr(*, "FUN")=function (x)
## ..- attr(*, "source")= chr "function (x) max(x, na.rm = TRUE)"
## - attr(*, "order.groups")= logi TRUE#Nhận xét sơ lược về khác biệt giữa nam và nữ
library(ggplot2)
theme_set(theme_bw())
ggplot(data = Orthodont) +
geom_point(aes(x = age, y = distance), col = "blue", size = 3) +
geom_line(aes(x = age, y = distance, group = Subject)) +
facet_grid(~ Sex)
#Vẽ đường thể hiện tương quan HQTT
lm = lm(distance ~ Sex:age + Sex, data = Orthodont)
summary(lm)##
## Call:
## lm(formula = distance ~ Sex:age + Sex, data = Orthodont)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.6156 -1.3219 -0.1682 1.3299 5.2469
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.3406 1.4162 11.538 < 2e-16 ***
## SexFemale 1.0321 2.2188 0.465 0.64279
## SexMale:age 0.7844 0.1262 6.217 1.07e-08 ***
## SexFemale:age 0.4795 0.1522 3.152 0.00212 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.257 on 104 degrees of freedom
## Multiple R-squared: 0.4227, Adjusted R-squared: 0.4061
## F-statistic: 25.39 on 3 and 104 DF, p-value: 2.108e-12Orthodont$pred.lm = predict(lm)
ggplot(data=Orthodont, add = T) +
geom_point(aes(x=age,y=distance), color="blue", size=3)+
geom_line(aes(x=age,y=pred.lm)) +
facet_grid(~Sex) 
2.2 Xây dựng mô hình đánh giá khoảng cách distance cho đối tượng i bằng những thông tin sau:
* yij: distance của đối tượng i ở tuổi j
* A0: trung bình distance của thiếu niên nam lúc 8 tuổi
* A1: khác biệt về distance giữa nam và nữ lúc 8 tuổi
* B0: tốc độ phát triển trung bình của distance ở thiếu niên nam
* B1: khác biệt về tốc độ phát triển của distance giữa nam và nữ
* u: dao động về distance giữa các đối tượng tham gia nghiên cứu lúc 8 tuổi
* v: dao độ tốc độ phát tiển distance giữa các đối tượng tham gia nghiên cứu
* eij: giá trị nhiễu (≠ giá trị đo được và giá trị dự báo từ mô hình) của người i.
- Sex (giới tính, 1= Nữ, 0= Nam);
- age (tuổi j ~ thời điểm đo lường distance)
Mô hình Unconditional Model : xét ảnh hưởng của age, không xét ảnh hưởng của Sex
Yij = Alphai + Betai + eij
* Alphai = A0 + random u
* Betai = B0 + random v
* Yij = (A0 + u) + (B0 + v).ageij + eij
* Yij = (A0 + B0.ageij) + (u + v.ageij + eij)
Mô hình Conditional Model : xét ảnh hưởng của age, xét thêm ảnh hưởng của Sex
* Alphai = A0 + A1.Sex + ui
* Betai = B0 + B1.Sex + vi
* Yij = (A0 + A1.Sex + u) + (B0 + B1.Sex + v).ageij + eij
* Yij = (A0 + A1.Sex + B0.age + B1.Sex.age) + (u + v.ageij + eij)
2.3 Thực hiện mô hình phân tích ảnh hưởng hỗn hợp bằng R
library(lme4)
m = lmer(distance ~ age + (1 + age|Subject), data = Orthodont)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.110724
## (tol = 0.002, component 1)summary(m)## Linear mixed model fit by REML ['lmerMod']
## Formula: distance ~ age + (1 + age | Subject)
## Data: Orthodont
##
## REML criterion at convergence: 442.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1252 -0.4933 0.0078 0.4577 3.9490
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 6.33344 2.5166
## age 0.05692 0.2386 -0.65
## Residual 1.68315 1.2974
## Number of obs: 108, groups: Subject, 27
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 16.76111 0.79203 21.162
## age 0.66019 0.07229 9.133
##
## Correlation of Fixed Effects:
## (Intr)
## age -0.853
## convergence code: 0
## Model failed to converge with max|grad| = 0.110724 (tol = 0.002, component 1)m1 = lmer(distance ~ Sex + age + Sex*age + (1 + age|Subject), data = Orthodont)
summary(m1)## Linear mixed model fit by REML ['lmerMod']
## Formula: distance ~ Sex + age + Sex * age + (1 + age | Subject)
## Data: Orthodont
##
## REML criterion at convergence: 432.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1694 -0.3862 0.0069 0.4454 3.8491
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 5.77441 2.4030
## age 0.03245 0.1801 -0.67
## Residual 1.71661 1.3102
## Number of obs: 108, groups: Subject, 27
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 16.34063 1.01824 16.048
## SexFemale 1.03210 1.59528 0.647
## age 0.78437 0.08598 9.123
## SexFemale:age -0.30483 0.13471 -2.263
##
## Correlation of Fixed Effects:
## (Intr) SexFml age
## SexFemale -0.638
## age -0.880 0.562
## SexFemale:g 0.562 -0.880 -0.6382.4 Viết lại mô hình 2.2 bằng các ước số từ kết quả mô hình
distance = 16.34 + 1.03.Sex + 0.78.age - 0.30.Sex.age
- Nhóm nam : 16.34 - 0.78.age
- Nhóm nữ: 17.07 - 1.08.age
+ A0 = 16.34 : trung bình distance của thiếu niên nam lúc 8 tuổi
+ A1 = 1.03 : khác biệt về distance giữa nam và nữ lúc 8 tuổi
+ B0 = 0.78 : tốc độ phát triển trung bình của distance ở thiếu niên nam
+ B1 = -0.3 : khác biệt về tốc độ phát triển của distance giữa nam và nữ
+ u = 5.77 : dao động về distance giữa các đối tượng tham gia nghiên cứu lúc 8 tuổi
+ v = 0.03 : dao độ tốc độ phát tiển distance giữa các đối tượng tham gia nghiên cứu
+ eij = 1.72 : giá trị nhiễu (≠ giá trị đo được và giá trị dự báo từ mô hình) của người i.
2.5 Bạn kết luận như thế nào về kết quả nghiên cứu?