R: 線形混合モデルの結果を多重比較 [統計]
線形混合モデルの結果を多重比較してみるメモ。 足軽日記さんのmultcomp研究(1) をほぼなぞったような内容である。
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)
コメント 0