Triển khai qua lme4

Tạo dữ liệu – dataframe

Score = c(3, 13, 13, 8, 11, 9, 12, 7, 16, 15, 18, 6, 21, 34, 26, 11, 24, 14, 21, 5, 17, 17, 23, 19, 7)
Group = c(rep("1", 10), rep("2", 6), rep("3", 9))
dat = data.frame(Group, Score)
dat
##    Group Score
## 1      1     3
## 2      1    13
## 3      1    13
## 4      1     8
## 5      1    11
## 6      1     9
## 7      1    12
## 8      1     7
## 9      1    16
## 10     1    15
## 11     2    18
## 12     2     6
## 13     2    21
## 14     2    34
## 15     2    26
## 16     2    11
## 17     3    24
## 18     3    14
## 19     3    21
## 20     3     5
## 21     3    17
## 22     3    17
## 23     3    23
## 24     3    19
## 25     3     7

Phân tích ANOVA truyền thống

anova = aov(Score ~ Group, data=dat)
summary(anova)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## Group        2  312.6  156.28   3.413 0.0512 .
## Residuals   22 1007.4   45.79                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Triển khai mô hình bằng hàm lmer

library(lme4)
## Loading required package: Matrix
fit = lmer(Score ~ 1 + (1 | Group), data=dat, REML=0)
summary(fit)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Score ~ 1 + (1 | Group)
##    Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    175.0    178.7    -84.5    169.0       22 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.64958 -0.65718  0.08166  0.52496  2.48793 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Group    (Intercept)  7.083   2.661   
##  Residual             45.797   6.767   
## Number of obs: 25, groups:  Group, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   15.150      2.059   7.358

No P value!

Tìm P value bằng Package khác

library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
fit = lmer(Score ~ 1 + (1 | Group), data=dat, REML=0) 
summary(fit)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Score ~ 1 + (1 | Group)
##    Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    175.0    178.7    -84.5    169.0       22 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.64958 -0.65718  0.08166  0.52496  2.48793 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Group    (Intercept)  7.083   2.661   
##  Residual             45.797   6.767   
## Number of obs: 25, groups:  Group, 3
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)   
## (Intercept)   15.150      2.059  2.942   7.358  0.00554 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mô hình là: Score = 15.15 + error (Variance = 45.8)

Triển khai mô hình mới (ảnh hưởng của Group)

fit = lmer(Score ~ 1 + Group + (1 | Group), data=dat, REML=0)
## boundary (singular) fit: see ?isSingular
summary(fit)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Score ~ 1 + Group + (1 | Group)
##    Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    173.4    179.4    -81.7    163.4       20 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1004 -0.4253  0.1050  0.6774  2.3104 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Group    (Intercept)  0.0     0.000   
##  Residual             40.3     6.348   
## Number of obs: 25, groups:  Group, 3
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   10.700      2.007 25.000   5.330 1.59e-05 ***
## Group2         8.633      3.278 25.000   2.634   0.0143 *  
## Group3         5.633      2.917 25.000   1.931   0.0648 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) Group2
## Group2 -0.612       
## Group3 -0.688  0.421
## convergence code: 0
## boundary (singular) fit: see ?isSingular

```

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.