[丢掉HLM软件]基于R的多层模型分析教程贴

bfzldh 发表于 2016-12-6 09:59
上面的教程已经说得很清楚了，你看看吧

好的，谢谢。

是否搞过分层分位回归？

jiandong4388 发表于 2016-12-11 08:18
是否搞过分层分位回归？

没有哎，自己用不到的都没有花时间学。我们专业喜欢用连续变量，即使是分类的，也要想法设法把它变成连续的。

感觉高大上，但是R确实有难度，值得挑战

齐威王王 发表于 2016-12-13 13:19
感觉高大上，但是R确实有难度，值得挑战

单纯只是为了分析数据的话，其实难度没那么大，先把数据输入、输出学会就入门了。

7.使用lme4程序包lmer函数拟合随机截距模型
#lme4函数界定模型的语句更简洁，对复杂模型的界定方法也更灵活。我们使用lme4程序包中的lmer函数拟合Model3.1，

1. Model3.7 <- lmer(geread~gevocab +(1|school), data = Achieve)

复制代码

#不同于lme函数，lmer函数中随机效应是用括号表示出来的，而不是专门用random来声明。上述语句的计算结果如下#6%lme4 1.0以后的版本发生了变化，“[gn]lmer now produces objects of class merMod rather than class mer as before”
#6%仅输入命令Model3.7，可得到以下结果

1. Model3.7

复制代码

Linear mixed model fit by REML ['lmerMod']
Formula: geread ~ gevocab + (1 | school)
   Data: Achieve
REML criterion at convergence: 43137.2
Random effects:
Groups      Name           Std.Dev.
school       (Intercept)    0.3159  
Residual                         1.9407  
Number of obs: 10320, groups:  school, 160
Fixed Effects:
(Intercept)      gevocab  
     2.0234       0.5129  
#6%输入summary(Model3.7),可得到以下结果：

1. summary(Model3.7)

复制代码

Linear mixed model fit by REML ['merModLmerTest']
Formula: geread ~ gevocab + (1 | school)
   Data: Achieve

REML criterion at convergence: 43137.2

Scaled residuals:
    Min      1Q          Median      3Q         Max
-3.0823   -0.5735  -0.2103     0.3207  4.4334

Random effects:
Groups     Name          Variance   Std.Dev.
school     (Intercept)   0.09978   0.3159  
Residual                     3.76647   1.9407  
Number of obs: 10320, groups:  school, 160

Fixed effects:
                     Estimate    Std. Error      t value
(Intercept)    2.023356   0.049309     41.03
gevocab        0.512898   0.008373     61.26

Correlation of Fixed Effects:
               (Intr)
gevocab -0.758
#6%可见，上述信息中不包含固定效应的显著性检验的相关信息（p值等）。
#尽管nlme同时报告了t值及其p值，但lme4的固定效应部分仅报告了t值。
#在lme中根据t分布确定p值这一方法并不能得到正确的p值。因此，有必要采用其他的方法确定p值。
#为了得到关于显著性的结果，可以使用lmerTest程序包中的summary函数。
#该函数基于Satterthwaite approximations对固定效应的显著性进行检验。

1. library(lmerTest)
   
2. summary(Model3.7)

复制代码

Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees of freedom [lmerMod]
Formula: geread ~ gevocab + (1 | school)
   Data: Achieve

REML criterion at convergence: 43137.2

Scaled residuals:
    Min      1Q       Median      3Q          Max
-3.0823 -0.5735 -0.2103      0.3207  4.4334

Random effects:
Groups   Name        Variance   Std.Dev.
school   (Intercept)  0.09978  0.3159  
Residual                  3.76647  1.9407  
Number of obs: 10320, groups:  school, 160

Fixed effects:
                  Estimate      Std. Error        df                  t value  Pr(>|t|)   
(Intercept) 2.023e+00  4.931e-02       7.580e+02    41.03   <2e-16 ***
gevocab     5.129e-01   8.373e-03        9.801e+03   61.26   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
              (Intr)
gevocab -0.758

8.使用lme4拟合零模型
#6%使用lme4拟合零模型。

1. Model3.7.0 <- lmer(geread~(1|school), data = Achieve)
   
2. summary(Model3.7.0)

复制代码

#结果如下
Linear mixed model fit by REML ['lmerMod']
Formula: geread ~ (1 | school)
   Data: Achieve

REML criterion at convergence: 46268.3

Scaled residuals:
    Min      1Q       Median      3Q        Max
-2.3229 -0.6378 -0.2138     0.2850  3.8812

Random effects:
Groups   Name         Variance    Std.Dev.
school   (Intercept)   0.3915     0.6257  
Residual                   5.0450     2.2461  
Number of obs: 10320, groups:  school, 160

Fixed effects:
                  Estimate    Std. Error   t value
(Intercept)  4.30675    0.05498    78.34

#使用lmerTest中的summary()函数

1. library(lmerTest)

复制代码

#结果显示，计算出错
summary from lme4 is returned
some computational error has occurred in lmerTest
……

9.使用lme4（lmerTest）在模型中加入多个自变量，或者在模型中加入交互作用#下面的例子包含了多个自变量：

1. library(lmerTest)
   
2. Model3.8 <- lmer(geread~gevocab + senroll +(1|school), data = Achieve)
   
3. summary(Model3.8)

复制代码

#结果如下

Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees of freedom [lmerMod]

Formula: geread ~ gevocab + senroll + (1 | school)

   Data: Achieve

REML criterion at convergence: 43152.1

Scaled residuals:

    Min      1Q       Median      3Q       Max

-3.0834 -0.5729  -0.2103     0.3212  4.4336

Random effects:

Groups   Name         Variance  Std.Dev.

school   (Intercept)  0.1003    0.3168  

Residual                  3.7665    1.9408  

Number of obs: 10320, groups:  school, 160

Fixed effects:

                   Estimate      Std. Error    df                 t value   Pr(>|t|)

(Intercept)  2.075e+00  1.140e-01   2.370e+02   18.20    <2e-16 ***

gevocab      5.129e-01   8.373e-03   9.798e+03   61.25    <2e-16 ***

senroll         -1.026e-04  2.051e-04   1.650e+02   -0.50     0.618   

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:

              (Intr)    gevocb

gevocab -0.327      

senroll   -0.901  -0.002

#下面在lme4(lmerTest)中加入交互项，这里加入水平1的gevocab和age的交互项，跨水平的交互作用请自行尝试

#如果不进行中心化，计算可能会失败

1. cgevocab=Achieve$gevocab-mean(Achieve$gevocab)
   
2. cage=Achieve$age-mean(Achieve$age)
   
3. Model3.8.s <-
   
4.   lmer(geread ~ cgevocab + cage + cgevocab*cage + (1 | school),
   
5.        data = Achieve)
   
6. summary(Model3.8.s)

复制代码

#结果如下：

Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees   of freedom [lmerMod]

Formula: geread ~ cgevocab + cage + cgevocab * cage + (1 | school)

   Data: Achieve

REML criterion at convergence: 43143.5

Scaled residuals:

    Min      1Q       Median      3Q         Max

-3.0635 -0.5706 -0.2108     0.3191  4.4467

Random effects:

Groups   Name        Variance   Std.Dev.

school   (Intercept) 0.09875   0.3143  

Residual                 3.76247    1.9397  

Number of obs: 10320, groups:  school, 160

Fixed effects:

                          Estimate      Std. Error     df                 t value     Pr(>|t|)   

(Intercept)         4.332e+00  3.206e-02    1.540e+02   135.124  < 2e-16 ***

cgevocab           5.125e-01   8.380e-03     9.819e+03   61.159    < 2e-16 ***

cage                  -6.777e-03  3.917e-03     1.032e+04  -1.730       0.08366 .  

cgevocab:cage  5.027e-03    1.750e-03     1.030e+04   2.873       0.00407 **

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:

                     (Intr)    cgevcb  cage

cgevocab       0.008              

cage              0.007   0.053      

cgevocab:cg  0.043   0.021   0.205

10.使用lme4(lmerTest)拟合随机斜率模型

#下面的模型仅包含一个水平1的自变量gevocab，并设定其斜率在水平2自由变化，只需在括号中声明gevocab

#只有提前中心化，才能成功计算

1. library(lmerTest)

复制代码

#结果如下

Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees of freedom [lmerMod]

Formula: geread ~ cgevocab + (cgevocab | school)

   Data: Achieve

REML criterion at convergence: 42992.9

Scaled residuals:

Min        1Q           Median   3Q         Max

-3.7102  -0.5674  -0.2074    0.3176  4.6775

Random effects:

Groups   Name          Variance Std.Dev. Corr

school   (Intercept)    0.1025    0.3202      

                cgevocab    0.0193    0.1389    0.52

Residual                    3.6659    1.9147      

Number of obs: 10320, groups:  school, 160

Fixed effects:

                   Estimate   Std. Error   df               t value  Pr(>|t|)   

(Intercept)  4.34411    0.03271     145.11000  132.8   <2e-16 ***

cgevocab    0.52036    0.01442     144.87000  36.1    <2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:

         (Intr)

cgevocab 0.362

11.在lme4中将多个随机效应设定为相关的或相互独立的
#在lme4中，我们可以界定多个自变量的随机斜率之间之间是否相关。
#下面呈现了Model3.10和Model3.11的语句，
#Model3.10的随机效应写在两个单独的括号里，这表示相应的两个随机效应是不相关的，
#Model3.11中，随机效应写在同一个括号里，表示相应的两个随机效应是相关的

1. Model3.10 <-
   
2. lmer(geread~cgevocab + cage+(cgevocab + cage|school), Achieve)
   
3. Model3.11 <-
   
4. lmer(geread~cgevocab + cage+ (cgevocab|school) + (cage|school),Achieve)

复制代码

1. summary(Model3.10)

复制代码

#Model3.10的结果如下
Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees of freedom [lmerMod]
Formula: geread ~ cgevocab + cage + (cgevocab + cage | school)
   Data: Achieve

REML criterion at convergence: 42995.3

Scaled residuals:
    Min      1Q      Median      3Q        Max
-3.6735 -0.5682 -0.2091    0.3184  4.6840

Random effects:
Groups   Name          Variance      Std.Dev.    Corr      
school   (Intercept)   1.022e-01   0.319698            
              cgevocab     1.902e-02   0.137918   0.53      
              cage            2.509e-05   0.005009  -0.28  -0.96
Residual                   3.664e+00   1.914181            
Number of obs: 10320, groups:  school, 160

Fixed effects:
                   Estimate      Std. Error         df                 t value     Pr(>|t|)   
(Intercept)  4.344e+00  3.269e-02         1.442e+02  132.896   <2e-16 ***
cgevocab     5.193e-01  1.435e-02          1.458e+02  36.187     <2e-16 ***
cage           -8.882e-03  3.822e-03          2.632e+03  -2.324       0.0202 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
                (Intr)   cgevcb
cgevocab  0.367      
cage        -0.020  -0.048
#Model3.11的结果如下

1. summary(Model3.11)

复制代码

Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees of freedom [lmerMod]
Formula: geread ~ cgevocab + cage + (cgevocab | school) + (cage | school)
   Data: Achieve

REML criterion at convergence: 42996.5

Scaled residuals:
    Min      1Q       Median      3Q        Max
-3.6937 -0.5681 -0.2081     0.3182  4.6744

Random effects:
Groups   Name         Variance       Std.Dev.       Corr
school    (Intercept)  2.914e-02   0.1707126     
               cgevocab    1.919e-02   0.1385421   1.00
school.1 (Intercept)  7.272e-02   0.2696680     
                cage          7.523e-07   0.0008674    1.00
Residual                    3.665e+00 1.9143975     
Number of obs: 10320, groups:  school, 160

Fixed effects:
                   Estimate      Std. Error         df                 t value      Pr(>|t|)   
(Intercept)  4.344e+00  3.264e-02         1.430e+02  133.065   <2e-16 ***
cgevocab     5.192e-01  1.440e-02          1.460e+02  36.063     <2e-16 ***
cage            -8.956e-03  3.798e-03         8.039e+03  -2.358      0.0184 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
                (Intr) cgevcb
cgevocab 0.368        
cage         0.015  0.033