LMMs
Nguyen Ngoc Thieu
October 15, 2018
Giới thiệu Linear Mixed Effect Models
Tôi muốn thực hiện một phân tích ANOVA nhưng data của tôi không thỏa mãn những điều kiện của ANOVA test. Tôi phải làm gì đây? Hoặc mô hình linear của tôi đã giải thích được 72.3% sự biến thiên của y rồi, tôi muốn thử tìm thêm biến số nào sẽ góp phần giải thích cho 27.7% còn lại, tôi phải làm gì nữa?
Dữ liệu tương quan rất thường gặp trong các phân tích thống kê. Điều này có thể là do việc phân nhóm các đối tượng nghiên cứu, ví dụ: những bệnh nhân trong các nhóm điều trị khác nhau, các nhóm tuổi khác nhau hoặc các phép đo lặp lại của từng đối tượng theo thời gian T0, T1, T2. Dữ liệu tạo thành một hệ thống phân cấp của các phép đo có tương quan với nhau (hierachy for levels of correlated measurements).
Phân tích mixed models cho ta một cách tiếp cận chung, rõ ràng trong các tình huống này, bởi vì nó cho phép một loạt các mô hình tương quan (hoặc cấu trúc variance-covariance) được mô hình hóa một cách rõ ràng.
Mixed models = fixed effect + random effect
Trong bộ dữ liệu dưới đây, chúng ta đo lường giá trị xét nghiệm y ở những bệnh nhân, từ 1-6, tại ba thời điểm T0, T1, T2. Những bệnh nhân này được điều trị bởi các thầy thuốc khác nhau (therapists: 1 - 4), theo hai loại điều trị A và B.
Sơ đồ phân cấp ba level theo thời gian như sau:
Fixed effect tác động lên trung bình biến số y trong các nhóm đối tượng
Random effect tác động lên variance của y trong các nhóm đối tượng
Nói theo ví dụ trên, biến số Treatment làm thay đổi y, tạo ra khác biệt của trung bình y giữa hai nhóm của Treatment.
Biến số Therapist (người điều trị) không phản ánh khác biệt có ý nghĩa của trung bình y trong các nhóm của Therapist, nhưng variance của y thì khác nhau. Sự khác biệt variance này cho thấy y biến động khác thường, có ý nghĩa thực tế, khi được điều trị bởi những Therapists khác nhau, mặc dù giá trị trung bình không khác nhau.
Mixed models cho ta cơ hội đánh giá đồng thời của hai loại tác động fixed và random lên y, điều này khác với ANOVA, chỉ đánh giá những biến số có fixed effects.
ANOVA vs. Mixed Effect Models
ANOVA: Đòi hỏi tất cả những đối tượng phải có đầy đủ dữ liệu cho số lần đo lường bằng nhau. VD mọi đối tượng phải có đủ 3 lần đo như nhau. Những đối tượng bị missing data hoặc có nhiều hơn, ít hơn 3 lần đo đều bị loại khỏi phân tích.
Những điều kiện (assumptions) của ANOVA:
- equal variances
- independence of errors
- normal distribution of errors
- additivity of treatment effects
Trong thực hành, những yêu cầu này thường khó đạt được hoàn toàn. Khi dữ liệu tương quan thì các điều kiện đối với ANOVA không còn thỏa mãn nữa. Khi các những điều kiện bị vi phạm thì kết quả phân tích ANOVA không còn giá trị nữa.
LMMs chấp nhận missing data và xem chúng như những trường hợp mất dữ liệu ngẫu nhiên và cho phép các đối tượng có số lần đo lường không bằng nhau (unbalanced clustered data).
Nhận biết một số categorical hoặc grouping variables là fixed hay random effects
Việc nhân biết và sắp xếp các biến số này vào mô hình với vai trò là fixed hay random effects tùy thuộc vào kinh nghiệm xử lí số liệu của người phân tích. Một số biến số thường gặp như sau:
- Fixed effect:
Nhóm điều trị, Giới tính, Centers, Khu vực, Địa phương
- Random effect:
Phân nhóm tuổi, Cá nhân đối tượng, Thời điểm theo dõi, Hộ gia đình, Lứa gia súc
Kiểm tra việc nhận diện, phân loại một biến số là fixed hoặc random đúng hay chưa bằng cách đánh giá mô hình LMMs để biết ý nghĩa thống kê của mô hình khi có sự tham gia của biến số đó.
Thực hành
Chuẩn bị dữ liệu
library(tidyverse)## -- Attaching packages ---------------------------------- tidyverse 1.2.1 --## v ggplot2 2.2.1 v purrr 0.2.4
## v tibble 1.4.1 v dplyr 0.7.4
## v tidyr 0.7.2 v stringr 1.2.0
## v readr 1.1.1 v forcats 0.2.0## -- Conflicts ------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(data.table)##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, last## The following object is masked from 'package:purrr':
##
## transposemydat <- fread('http://www-personal.umich.edu/~bwest/rat_pup.dat')
head(mydat,15)## pup_id weight sex litter litsize treatment
## 1: 1 6.60 Male 1 12 Control
## 2: 2 7.40 Male 1 12 Control
## 3: 3 7.15 Male 1 12 Control
## 4: 4 7.24 Male 1 12 Control
## 5: 5 7.10 Male 1 12 Control
## 6: 6 6.04 Male 1 12 Control
## 7: 7 6.98 Male 1 12 Control
## 8: 8 7.05 Male 1 12 Control
## 9: 9 6.95 Female 1 12 Control
## 10: 10 6.29 Female 1 12 Control
## 11: 11 6.77 Female 1 12 Control
## 12: 12 6.57 Female 1 12 Control
## 13: 13 6.37 Male 2 14 Control
## 14: 14 6.37 Male 2 14 Control
## 15: 15 6.90 Male 2 14 Controldim(mydat)## [1] 322 6Bộ dữ liệu về weight của những con chuột con, dataset có 322 dòng và 6 biến số như trên.Trong đó treatment là biến số nhóm điều trị (high, low và control), litter là bầy chuột được sinh ra bởi một chuột mẹ, litsize là số con trong một bầy chuột nhất định.
Mục đích của nghiên cứu là tìm quy luật biến đổi của trọng lượng chuột con (weight) theo giới tính, số con trong một lứa và dưới tác động của loại điều trị.
weight ~ treatment + sex + litter + treatment:sex + treatment:litter
Tóm tắt về weight của chuột con trong các nhóm
library(plyr)## -------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## -------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following object is masked from 'package:purrr':
##
## compactddply(mydat,~treatment*sex,summarise,mean=mean(weight),sd=sd(weight),min=min(weight),max=max(weight))## treatment sex mean sd min max
## 1 Control Female 6.116111 0.6851179 3.68 7.57
## 2 Control Male 6.471039 0.7537880 4.57 8.33
## 3 High Female 5.851562 0.6001887 4.48 7.68
## 4 High Male 5.918485 0.6909058 5.01 7.70
## 5 Low Female 5.837538 0.4504964 4.75 7.73
## 6 Low Male 6.025082 0.3803403 5.25 7.13table(mydat$treatment,mydat$sex)##
## Female Male
## Control 54 77
## High 32 33
## Low 65 61p<-ggplot(mydat, aes(x=treatment, y=weight, fill=sex)) +
geom_boxplot(position=position_dodge(1))+
ggtitle("Weight By Treatment Groups and Sex")+
stat_summary(fun.y="mean", geom="point", shape=23, size=4)
p
Nhìn vào bảng tóm tắt và biểu đồ, ta có những nhận xét sơ bộ như sau: trọng lượng (weight) của chuột con càng thấp khi có can thiệp (high treatment vs low treatment) so với không can thiệp (control) và chuột đực nặng hơn chuột cái.
Ngoài ra, nhìn vào biểu đồ có hình thoi (mean) và box (median) ta thấy variance của weight trong các nhóm tương đối khác nhau, đặt ra nghi vấn về unequal variances.
ANOVA: kiểm tra các điều kiện (assumptions)
Trước hết chúng ta phân tích theo mô hình ANOVA xem weight thay đổi theo các nhóm điều trị, giới tính, số con trong lứa sinh như thế nào.
# Chuẩn bị dữ liệu cho ANOVA test
ratpup <- mydat
ratpup$sex1[ratpup$sex == "Female"] <- 1
ratpup$sex1[ratpup$sex == "Male"] <- 0
head(ratpup)## pup_id weight sex litter litsize treatment sex1
## 1: 1 6.60 Male 1 12 Control 0
## 2: 2 7.40 Male 1 12 Control 0
## 3: 3 7.15 Male 1 12 Control 0
## 4: 4 7.24 Male 1 12 Control 0
## 5: 5 7.10 Male 1 12 Control 0
## 6: 6 6.04 Male 1 12 Control 0# Anova model
fit <- aov(weight ~ treatment + sex1 + litter + treatment:sex1, data=ratpup)
summary(fit)## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 2 13.20 6.601 17.927 4.23e-08 ***
## sex1 1 4.19 4.194 11.389 0.000831 ***
## litter 1 0.16 0.156 0.424 0.515301
## treatment:sex1 2 1.01 0.504 1.368 0.256015
## Residuals 315 115.99 0.368
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Weight liên quan có ý nghĩa với treatment và sex1, không liên quan có ý nghĩa với litter và tương tác treatment:sex1.
# kiểm tra tính normal distribution of errors
fit2 <- aov(weight ~ treatment + sex1, data=ratpup)
ratpup$resids <- residuals(fit2)
ratpup$preds <- predict(fit2)
ratpup$sq_preds <- ratpup$preds^2
plot(resids ~ preds, data = ratpup, xlab = "Predicted Values", ylab = "Residuals") 
shapiro.test(ratpup$resids) ##
## Shapiro-Wilk normality test
##
## data: ratpup$resids
## W = 0.97819, p-value = 8.285e-05p-value = 8.285e-05 cho thấy residuals (errors) không có phân phối normal.
#Perform Levene's Test for homogenity of variances
library(car) ## Warning: package 'car' was built under R version 3.4.4## Loading required package: carData## Warning: package 'carData' was built under R version 3.4.4##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## someleveneTest(weight ~ treatment, data = ratpup) ## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 10.921 2.586e-05 ***
## 319
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(weight ~ sex, data = ratpup)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.2755 0.07126 .
## 320
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1p = 2.586e-05 *** cho thấy variance trong các nhóm của treatment không đồng nhất (unequal variance). Tôi sẽ trở lại vấn đề này trong phần sau.
Như vậy có ít nhất 2 assumptions của ANOVA test bị vi phạm, cho nên kết quả ANOVA test sẽ không có giá trị.
Việc kiểm tra những assumptions này có thể bị bỏ qua trong thực hành phân tích dữ liệu, dẫn đến kết quả phân tích không đáng tin cậy.
Đó là lí do chúng ta cần phải xem xét đến mixed effect models cho những trường hợp như thế này.
Linear Mixed Effect Models với lme()
library(nlme)## Warning: package 'nlme' was built under R version 3.4.4##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapsemodel3.1.fit <- lme(weight ~ treatment + sex1 + litsize +
treatment:sex1, random = ~ 1 | litter,
data = ratpup, method = "REML")
summary(model3.1.fit)## Linear mixed-effects model fit by REML
## Data: ratpup
## AIC BIC logLik
## 419.1043 452.8775 -200.5522
##
## Random effects:
## Formula: ~1 | litter
## (Intercept) Residual
## StdDev: 0.3106722 0.404337
##
## Fixed effects: weight ~ treatment + sex1 + litsize + treatment:sex1
## Value Std.Error DF t-value p-value
## (Intercept) 8.323340 0.27333009 292 30.451605 0.0000
## treatmentHigh -0.906057 0.19154238 23 -4.730320 0.0001
## treatmentLow -0.467040 0.15818328 23 -2.952521 0.0071
## sex1 -0.411688 0.07315410 292 -5.627679 0.0000
## litsize -0.128382 0.01875336 23 -6.845819 0.0000
## treatmentHigh:sex1 0.107023 0.13176318 292 0.812239 0.4173
## treatmentLow:sex1 0.083866 0.10568189 292 0.793568 0.4281
## Correlation:
## (Intr) trtmnH trtmnL sex1 litsiz trtH:1
## treatmentHigh -0.562
## treatmentLow -0.297 0.404
## sex1 -0.111 0.158 0.191
## litsize -0.916 0.363 0.022 0.001
## treatmentHigh:sex1 0.043 -0.323 -0.106 -0.555 0.020
## treatmentLow:sex1 0.044 -0.096 -0.320 -0.692 0.036 0.385
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -7.47250744 -0.50014749 0.02911668 0.57348178 3.00962055
##
## Number of Observations: 322
## Number of Groups: 27Phần fixed effect: treatment + sex1 + litsize + treatment:sex1
Phần random effect: random = ~ 1 | litter
Giải thích kết quả:
Các chỉ số AIC, BIC, LogLik nói lên performance của model này. Chúng ta dùng các chỉ số này để so sánh với các models biến đổi sau này nhằm tìm ra model tối ưu.
Intercept = 8.323340 là chỉ số weight trung bình.
Chuột trong nhóm “hight treatment” giảm trọng lượng -0.906057 đơn vị,
Chuột trong nhóm “low treatment” giảm trọng lượng -0.467040 đơn vị,
Chuột cái (sex = 1) nhẹ hơn chuột đực (sex = 0) -0.411688,
Trong những lứa (bầy) có số con tăng 1 thì trọng lượng trung bình của chuột giảm -0.128382 đơn vị.
Khi chọn biến số litter để đưa vào phần random effect trong mô hình, nó cho thấy có vai trò trong biến thiên của weight với SD = 0.3106722 và SD của residuals = 0.404337. Vai trò của biến số litter không được ANOVA test phát hiện vì thuật toán của ANOVA chỉ phát hiện fixed effects. Linear Mixed Effect Model (lme) giúp phát hiện ra sự khác biệt này.
Loại bỏ thành phần treatment:sex1 không có ý nghĩa thì mô hình còn lại sẽ là:
library(nlme)
model3.1b.fit <- lme(weight ~ treatment + sex1 + litsize , random = ~ 1 | litter,
data = ratpup, method = "REML")Package nlme cũng có hàm gls() cho ta phần fixed effect mà thôi. Chúng ta thử so sánh hai models có và không có phần random effect để xem random variable có ý nghĩa trong sự biến thiên cùa weight hay không:
# Model with fixed effects only
model3.lc.fit <- gls(weight ~ treatment + sex1 + litsize , data = ratpup)
anova(model3.1b.fit, model3.lc.fit)## Model df AIC BIC logLik Test L.Ratio p-value
## model3.1b.fit 1 7 410.9994 437.3117 -198.4997
## model3.lc.fit 2 6 499.4820 522.0354 -243.7410 1 vs 2 90.48262 <.0001Kết quả cho thấy model3.lc.fit có các chỉ số AIC, BIC, LogLik không ưu thế bằng mô hình model3.1b.fit (có phần random effect) và sự khác biệt giữa hai mô hình có ý nghĩa thống kê, p < 0.001.
Điều này nói lên vai trò có ý nghĩa của random effect (random = ~ 1 | litter) trong sự biến thiên của weight.
Đánh giá performance của model này qua biểu đồ residual và QQ như sau:
library(lattice)
plot(model3.1b.fit,resid(., type = "p") ~ fitted(.) | factor(treatment), layout = c(3,1),aspect=2, abline=0)
qqnorm(model3.1b.fit, ~ resid(.) | factor(treatment),layout = c(3,1) ,aspect = 2)
Variance của weight trong các nhóm high, low và control của treatment không đồng nhất, biểu lộ rõ trong biểu đồ QQ ở trên với khác biệt rõ trong nhóm control.
Còn nhớ ở phần trên, ta đã phát hiện các nhóm control, high và low treatment có variance không đồng nhất. Tôi sẽ làm rõ những thông số residual variance ở các levels khác nhau của biến số treatment như sau:
model3.2a.fit <- lme(weight ~ treatment + sex1 + litsize,
random = ~1 | litter, ratpup, method = "REML",
weights = varIdent(form = ~1 | treatment))
summary(model3.2a.fit)## Linear mixed-effects model fit by REML
## Data: ratpup
## AIC BIC logLik
## 373.1646 406.9947 -177.5823
##
## Random effects:
## Formula: ~1 | litter
## (Intercept) Residual
## StdDev: 0.3136624 0.5144278
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | treatment
## Parameter estimates:
## Control Low High
## 1.0000000 0.5643872 0.6354988
## Fixed effects: weight ~ treatment + sex1 + litsize
## Value Std.Error DF t-value p-value
## (Intercept) 8.322966 0.27310901 294 30.474886 0.0000
## treatmentHigh -0.862368 0.18297383 23 -4.713067 0.0001
## treatmentLow -0.433563 0.15163164 23 -2.859319 0.0089
## sex1 -0.343458 0.04180668 294 -8.215389 0.0000
## litsize -0.130331 0.01848141 23 -7.051998 0.0000
## Correlation:
## (Intr) trtmnH trtmnL sex1
## treatmentHigh -0.606
## treatmentLow -0.351 0.464
## sex1 -0.103 -0.005 -0.035
## litsize -0.913 0.399 0.066 0.044
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.97300195 -0.52157120 0.02077248 0.56022844 2.53765137
##
## Number of Observations: 322
## Number of Groups: 27Thành tố weights = varIdent(form = ~ 1 | treatment) khai báo cho model biết ta đang khai thác biệt variance trong các levels của treatment.
Trong phần random effects cho thấy Residual SD = 0.5144278
trong phần Parameter estimates:
Control | Low | High
1.0000000 | 0.5643872 | 0.6354988
Ta nhân Residual SD = 0.5144278 cho từng parameter estimates ở trên để cho ra Residual SD cho từng levels của treatment.
Như vậy SD cho mỗi nhóm treatment là:
Control: 0.5144
Low: 0.2901
high: 0.3266
Như vậy SD trong nhóm Control lớn nhất, nói cách khác, weight trong nhóm Control biến động nhiều hơn các nhóm khác.
So sánh hai model:
anova(model3.1b.fit, model3.2a.fit)## Model df AIC BIC logLik Test L.Ratio p-value
## model3.1b.fit 1 7 410.9994 437.3117 -198.4997
## model3.2a.fit 2 9 373.1646 406.9947 -177.5823 1 vs 2 41.83483 <.0001Vậy model có sự xem xét của treatment variance cho kết quả tốt hơn.
Ngoài ra chúng ta có thể thay đổi phương pháp từ method = “REML” sang method = “ML” để so sánh và tìm ra models tối ưu, khi so sánh các chỉ số AIC, BIC và LogLik của chúng.
Hàm lmer() từ package lme4
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following object is masked from 'package:tidyr':
##
## expand##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmListmodel3.1.fit.lmer <- lmer(weight ~ treatment + sex1 + litsize +
(1 | litter), ratpup, REML = TRUE)Trong hàm lmer() phần random có cách viết hơi khác với hàm lme() một tí, chúng ta chú ý.
summary(model3.1.fit.lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: weight ~ treatment + sex1 + litsize + (1 | litter)
## Data: ratpup
##
## REML criterion at convergence: 397
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.5743 -0.4699 0.0111 0.5712 3.0603
##
## Random effects:
## Groups Name Variance Std.Dev.
## litter (Intercept) 0.0974 0.3121
## Residual 0.1628 0.4035
## Number of obs: 322, groups: litter, 27
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8.30987 0.27371 30.360
## treatmentHigh -0.85870 0.18181 -4.723
## treatmentLow -0.42850 0.15040 -2.849
## sex1 -0.35908 0.04749 -7.562
## litsize -0.12900 0.01879 -6.864
##
## Correlation of Fixed Effects:
## (Intr) trtmnH trtmnL sex1
## treatmntHgh -0.581
## treatmentLw -0.300 0.423
## sex1 -0.111 -0.009 -0.040
## litsize -0.920 0.390 0.035 0.043Kết quả được trình bày tương tự hàm lme(). Tuy nhiên, không có p values.
Để xem phần random effect mà thôi chúng ta dùng hàm ranef()
ranef(model3.1.fit.lmer)## $litter
## (Intercept)
## 1 0.17805612
## 2 -0.07071525
## 3 -0.18289450
## 4 -0.05418319
## 5 0.34615302
## 6 -0.05363520
## 7 0.39046482
## 8 -0.02968245
## 9 -0.60785892
## 10 0.08429556
## 11 0.03832949
## 12 0.02490328
## 13 -0.36672606
## 14 0.01557233
## 15 0.03327425
## 16 0.07051165
## 17 -0.39760275
## 18 0.43692016
## 19 -0.18154659
## 20 0.32636424
## 21 -0.29062992
## 22 -0.49231419
## 23 0.26451251
## 24 0.18828098
## 25 0.23099260
## 26 0.24197633
## 27 -0.14281830Muốn thể hiện các giá trị p-values ta gọi package lmerTest trước khi viết model như sau:
library(lmerTest)## Warning: package 'lmerTest' was built under R version 3.4.4##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## stepmodel3.1.fit.lmer <- lmer(weight ~ treatment + sex1 + litsize +
(1 | litter), ratpup, REML = TRUE)
summary(model3.1.fit.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight ~ treatment + sex1 + litsize + (1 | litter)
## Data: ratpup
##
## REML criterion at convergence: 397
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.5743 -0.4699 0.0111 0.5712 3.0603
##
## Random effects:
## Groups Name Variance Std.Dev.
## litter (Intercept) 0.0974 0.3121
## Residual 0.1628 0.4035
## Number of obs: 322, groups: litter, 27
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.30987 0.27371 32.61000 30.360 < 2e-16 ***
## treatmentHigh -0.85870 0.18181 24.98000 -4.723 7.65e-05 ***
## treatmentLow -0.42850 0.15040 22.90000 -2.849 0.0091 **
## sex1 -0.35908 0.04749 301.82000 -7.562 4.81e-13 ***
## litsize -0.12900 0.01879 31.67000 -6.864 9.63e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmnH trtmnL sex1
## treatmntHgh -0.581
## treatmentLw -0.300 0.423
## sex1 -0.111 -0.009 -0.040
## litsize -0.920 0.390 0.035 0.043Hy vọng bài viết sẽ hữu ích với các bạn, giúp giải đáp những hạn chế của ANOVA test trong thực hành phân tích dữ liệu.
Cám ơn.
KÌ TỚI:
Linear Mixed Effect Models cho dữ liệu longitudinal, tức là phép đo lặp lại ở nhiều thời điểm. ______________________________________________________________