R: 線形混合モデルの結果を多重比較

線形混合モデルの結果を多重比較してみるメモ。 足軽日記さんのmultcomp研究（１） をほぼなぞったような内容である。

nlmeパッケージを使用。

library(nlme)

データを用意する。

set.seed(123)
n.block <- 10
n.samp <- 20
block <- factor(rep(1:n.block, n.samp))
treat <- factor(rep(c("a", "b", "c"), each = n.block * n.samp))
ranef <- rnorm(n.block, 0, 1)
y <- 10 + 3 * as.numeric(treat) + ranef[block] + rnorm(length(treat), 0, 2)
data <- data.frame(treat = treat, block = block, y = y)
pairs(data)

lme()にて解析。blockを切片に対するランダム効果に。

fit <- lme(y ~ treat, random = ~1 | block, data = data)
summary(fit)

## Linear mixed-effects model fit by REML
##  Data: data 
##    AIC  BIC logLik
##   2534 2556  -1262
## 
## Random effects:
##  Formula: ~1 | block
##         (Intercept) Residual
## StdDev:      0.9723    1.936
## 
## Fixed effects: y ~ treat 
##              Value Std.Error  DF t-value p-value
## (Intercept) 13.081    0.3366 588   38.87       0
## treatb       3.019    0.1936 588   15.59       0
## treatc       6.054    0.1936 588   31.27       0
##  Correlation: 
##        (Intr) treatb
## treatb -0.288       
## treatc -0.288  0.500
## 
## Standardized Within-Group Residuals:
##       Min        Q1       Med        Q3       Max 
## -3.176575 -0.619689 -0.008971  0.678863  3.356871 
## 
## Number of Observations: 600
## Number of Groups: 10

multcompパッケージのglht()関数にて多重比較。

library(multcomp)

## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: splines
## Loading required package: TH.data

summary(glht(fit, linfct = mcp(treat = "Tukey")))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lme.formula(fixed = y ~ treat, data = data, random = ~1 | block)
## 
## Linear Hypotheses:
##            Estimate Std. Error z value Pr(>|z|)    
## b - a == 0    3.019      0.194    15.6   <2e-16 ***
## c - a == 0    6.054      0.194    31.3   <2e-16 ***
## c - b == 0    3.035      0.194    15.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

タグ：多重比較 R