KU-R Session 3

本日の内容

検定

• 連続変数: t.test(), wilcox.test(), oneway.test(), kruskal.test()
• カテゴリカル変数: chisq.test(), fisher.test()
• 生存時間変数: survfit()とカプランマイヤー曲線, survdiff()

回帰分析

• 線形回帰モデル(Ordinary least square linear regression)
• 一般化線形回帰モデル(Generalized linear models)
• 生存分析
• 線形混合効果モデル(Linear mixed effects model)
• 一般化推定方程式(Generalized estimating equation)

グラフ作成

• 利用可能なグラフシステムの種類
• baseグラフィックによる作図
• ggplot2パッケージによる作図
• 3Dグラフ
• グラフの出力

検定

連続変数: t.test(), wilcox.test(), oneway.test(), kruskal.test()

## データ読み込み
nhefs <- XLConnect::readWorksheetFromFile("./nhefs_book.xls", sheet = 1)

## 二群
## 正規分布
## デフォルトでは両群の分散の一致を仮定しない検定になります。
t.test(age ~ sex, data = nhefs)

## 
##  Welch Two Sample t-test
## 
## data:  age by sex
## t = 1.7044, df = 1735.968, p-value = 0.08848
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1487584  2.1225396
## sample estimates:
## mean in group 0 mean in group 1 
##        44.13542        43.14853

## 通常のt検定はこちら
t.test(age ~ sex, data = nhefs, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  age by sex
## t = 1.7053, df = 1744, p-value = 0.08833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1481994  2.1219806
## sample estimates:
## mean in group 0 mean in group 1 
##        44.13542        43.14853

## 正規性のチェックはQQ plotが良いでしょうか。
car::qqPlot(nhefs$age)

## 非正規分布
## Wilcoxon test()
wilcox.test(age ~ sex, data = nhefs)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  age by sex
## W = 397971.5, p-value = 0.1075
## alternative hypothesis: true location shift is not equal to 0

## 多群
## ANOVAはoneway.test()で出来ます。
## これもデフォルトでは分散の一致を仮定しないものになります。
oneway.test(age ~ alcoholtype, data = nhefs)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  age and alcoholtype
## F = 32.5121, num df = 3.000, denom df = 579.648, p-value < 2.2e-16

oneway.test(age ~ alcoholtype, data = nhefs, var.equal = TRUE)

## 
##  One-way analysis of means
## 
## data:  age and alcoholtype
## F = 31.3619, num df = 3, denom df = 1742, p-value < 2.2e-16

## 非正規の場合はkruskal.testが使えます
kruskal.test(age ~ alcoholtype, data = nhefs)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  age by alcoholtype
## Kruskal-Wallis chi-squared = 95.596, df = 3, p-value < 2.2e-16

ANOVAに関しては実際は回帰分析の関数を使って行う方が一般的です。平方和のテーブルなどが出ます。

lm1 <- lm(age ~ factor(alcoholtype), data = nhefs)
anova(lm1)

## Analysis of Variance Table
## 
## Response: age
##                       Df Sum Sq Mean Sq F value    Pr(>F)    
## factor(alcoholtype)    3  13086  4361.9  31.362 < 2.2e-16 ***
## Residuals           1742 242284   139.1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Rでは結果がオブジェクトとして返ってきます。名前を付けて保存して内容をチェックします。

resTTest <- t.test(age ~ sex, data = nhefs)
## 名前だけ評価すると適切な方法で表示されます
resTTest

## 
##  Welch Two Sample t-test
## 
## data:  age by sex
## t = 1.7044, df = 1735.968, p-value = 0.08848
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1487584  2.1225396
## sample estimates:
## mean in group 0 mean in group 1 
##        44.13542        43.14853

## 要素の名前をチェックします。
names(resTTest)

## [1] "statistic"   "parameter"   "p.value"     "conf.int"    "estimate"   
## [6] "null.value"  "alternative" "method"      "data.name"

## p-valueだけ抜きだします。
resTTest$p.value

## [1] 0.08848213

カテゴリカル変数: chisq.test(), fisher.test()

カテゴリカル変数では前述の関数を使ってクロス集計表を作成して検定します。

## 表作成
tabSexQsmk <- xtabs( ~ sex + qsmk, data = nhefs)

## カイ二乗検定 デフォルトでは連続性補正をします。
chisq.test(tabSexQsmk)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabSexQsmk
## X-squared = 8.4915, df = 1, p-value = 0.003568

## 通常のカイ二乗検定
chisq.test(tabSexQsmk, correct = FALSE)

## 
##  Pearson's Chi-squared test
## 
## data:  tabSexQsmk
## X-squared = 8.8102, df = 1, p-value = 0.002996

## すでに表のデータの場合は下記のように入力しても良いでしょう。
## 4つの値を与えて、2列の行列を作成しています。標準では縦方向に数字が入ります。
## byrowで横向きに変更できます。
mat <- matrix(c(5,1,10,20), ncol = 2, byrow = TRUE)
resFisher <- fisher.test(mat)
resFisher

## 
##  Fisher's Exact Test for Count Data
## 
## data:  mat
## p-value = 0.06281
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##    0.8811891 494.4191744
## sample estimates:
## odds ratio 
##   9.364462

## 同じようにオブジェクトが返っています。
names(resFisher)

## [1] "p.value"     "conf.int"    "estimate"    "null.value"  "alternative"
## [6] "method"      "data.name"

## 2x2表のオッズ比を抜き出してみます。
resFisher$estimate

## odds ratio 
##   9.364462

生存時間変数: survfit()とカプランマイヤー曲線, survdiff()

デモデータが生存分析データではないので別のデータを使います。

## survivalパッケージから白血病のデータを読み込みます。
library(survival)
data(leukemia)
## チェック
head(leukemia)

##   time status          x
## 1    9      1 Maintained
## 2   13      1 Maintained
## 3   13      0 Maintained
## 4   18      1 Maintained
## 5   23      1 Maintained
## 6   28      0 Maintained

## カプランーマイヤーはこのように行います。
## 群分けしない場合はformulaの右には1とします (モデルに説明変数がなく切片のみという意味)
## 信頼区間のタイプはlog-log変換にしておいた方が狭くなります。
surv1 <- survfit(Surv(time, status) ~ 1, data = leukemia, conf.type = "log-log")

## 標準の結果表示では限られた情報のみ
surv1

## Call: survfit(formula = Surv(time, status) ~ 1, data = leukemia, conf.type = "log-log")
## 
## records   n.max n.start  events  median 0.95LCL 0.95UCL 
##      23      23      23      18      27      13      34

## summaryするとイベント時間ごとの集計がみられます。
summary(surv1)

## Call: survfit(formula = Surv(time, status) ~ 1, data = leukemia, conf.type = "log-log")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5     23       2   0.9130  0.0588      0.69495        0.978
##     8     21       2   0.8261  0.0790      0.60061        0.931
##     9     19       1   0.7826  0.0860      0.55421        0.903
##    12     18       1   0.7391  0.0916      0.50921        0.873
##    13     17       1   0.6957  0.0959      0.46564        0.842
##    18     14       1   0.6460  0.1011      0.41395        0.805
##    23     13       2   0.5466  0.1073      0.31925        0.726
##    27     11       1   0.4969  0.1084      0.27557        0.684
##    30      9       1   0.4417  0.1095      0.22738        0.637
##    31      8       1   0.3865  0.1089      0.18284        0.587
##    33      7       1   0.3313  0.1064      0.14183        0.535
##    34      6       1   0.2761  0.1020      0.10444        0.480
##    43      5       1   0.2208  0.0954      0.07100        0.422
##    45      4       1   0.1656  0.0860      0.04211        0.360
##    48      2       1   0.0828  0.0727      0.00696        0.287

## 任意の時間を指定するにはtimesを使います。
summary(surv1, times = c(0, 20, 40, 60))

## Call: survfit(formula = Surv(time, status) ~ 1, data = leukemia, conf.type = "log-log")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0     23       0   1.0000  0.0000      1.00000        1.000
##    20     13       8   0.6460  0.1011      0.41395        0.805
##    40      5       7   0.2761  0.1020      0.10444        0.480
##    60      1       3   0.0828  0.0727      0.00696        0.287

## KMはこの結果オブジェクトをplotすると描画されます。
plot(surv1)

## 信頼区間が不要な場合
plot(surv1, conf.int = FALSE)

## 群分けする場合はformulaを変えます。
surv2 <- survfit(Surv(time, status) ~ x, data = leukemia, conf.type = "log-log")
surv2

## Call: survfit(formula = Surv(time, status) ~ x, data = leukemia, conf.type = "log-log")
## 
##                 records n.max n.start events median 0.95LCL 0.95UCL
## x=Maintained         11    11      11      7     31      13      NA
## x=Nonmaintained      12    12      12     11     23       5      33

summary(surv2)

## Call: survfit(formula = Surv(time, status) ~ x, data = leukemia, conf.type = "log-log")
## 
##                 x=Maintained 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     9     11       1    0.909  0.0867       0.5081        0.987
##    13     10       1    0.818  0.1163       0.4474        0.951
##    18      8       1    0.716  0.1397       0.3502        0.899
##    23      7       1    0.614  0.1526       0.2658        0.835
##    31      5       1    0.491  0.1642       0.1673        0.753
##    34      4       1    0.368  0.1627       0.0928        0.657
##    48      2       1    0.184  0.1535       0.0117        0.525
## 
##                 x=Nonmaintained 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5     12       2   0.8333  0.1076      0.48171        0.956
##     8     10       2   0.6667  0.1361      0.33702        0.860
##    12      8       1   0.5833  0.1423      0.27014        0.801
##    23      6       1   0.4861  0.1481      0.19188        0.730
##    27      5       1   0.3889  0.1470      0.12627        0.650
##    30      4       1   0.2917  0.1387      0.07240        0.561
##    33      3       1   0.1944  0.1219      0.03120        0.461
##    43      2       1   0.0972  0.0919      0.00575        0.349
##    45      1       1   0.0000     NaN           NA           NA

summary(surv2, times = c(0, 20, 40, 60))

## Call: survfit(formula = Surv(time, status) ~ x, data = leukemia, conf.type = "log-log")
## 
##                 x=Maintained 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0     11       0    1.000   0.000       1.0000        1.000
##    20      7       3    0.716   0.140       0.3502        0.899
##    40      3       3    0.368   0.163       0.0928        0.657
##    60      1       1    0.184   0.153       0.0117        0.525
## 
##                 x=Nonmaintained 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0     12       0    1.000   0.000       1.0000        1.000
##    20      6       5    0.583   0.142       0.2701        0.801
##    40      2       4    0.194   0.122       0.0312        0.461

plot(surv2)

## log-rank検定の書式はsurvfitと同じです。
survdiff(Surv(time, status) ~ x, data = leukemia)

## Call:
## survdiff(formula = Surv(time, status) ~ x, data = leukemia)
## 
##                  N Observed Expected (O-E)^2/E (O-E)^2/V
## x=Maintained    11        7    10.69      1.27       3.4
## x=Nonmaintained 12       11     7.31      1.86       3.4
## 
##  Chisq= 3.4  on 1 degrees of freedom, p= 0.0653

## rmsパッケージのsurvplotの方がましなグラフが描けます。
library(rms)

## Loading required package: Hmisc
## Loading required package: grid
## Loading required package: lattice
## Loading required package: Formula
## Loading required package: ggplot2
## 
## Attaching package: 'Hmisc'
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, round.POSIXt, trunc.POSIXt, units
## 
## Loading required package: gridExtra
## Loading required package: SparseM
## 
## Attaching package: 'SparseM'
## 
## The following object is masked from 'package:base':
## 
##     backsolve

## npsurv()関数の書式はsurvfitと全く同じです。
np1 <- npsurv(Surv(time, status) ~ x, data = leukemia, conf.type = "log-log")
summary(np1)

## Call: npsurv(formula = Surv(time, status) ~ x, data = leukemia, conf.type = "log-log")
## 
##                 x=Maintained 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     9     11       1    0.909  0.0867       0.5081        0.987
##    13     10       1    0.818  0.1163       0.4474        0.951
##    18      8       1    0.716  0.1397       0.3502        0.899
##    23      7       1    0.614  0.1526       0.2658        0.835
##    31      5       1    0.491  0.1642       0.1673        0.753
##    34      4       1    0.368  0.1627       0.0928        0.657
##    48      2       1    0.184  0.1535       0.0117        0.525
## 
##                 x=Nonmaintained 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5     12       2   0.8333  0.1076      0.48171        0.956
##     8     10       2   0.6667  0.1361      0.33702        0.860
##    12      8       1   0.5833  0.1423      0.27014        0.801
##    23      6       1   0.4861  0.1481      0.19188        0.730
##    27      5       1   0.3889  0.1470      0.12627        0.650
##    30      4       1   0.2917  0.1387      0.07240        0.561
##    33      3       1   0.1944  0.1219      0.03120        0.461
##    43      2       1   0.0972  0.0919      0.00575        0.349
##    45      1       1   0.0000     NaN           NA           NA

## デフォルトのグラフ
survplot(np1)

## オプションの指定方法。
survplot(fit             = np1,
         conf            = c("none","bands","bars")[1],
         xlab            = "", ylab = "Survival",
         ## xlim         = c(0,100),
         label.curves    = TRUE,                  # label curves directly
         ## label.curves = list(keys = "lines"),  # legend instead of direct label
         levels.only     = TRUE,                  # show only levels, no variable name
         abbrev.label    = FALSE,                 # if label used, abbreviate
         ## fun          = function(x) {1 - x},   # Cumulative probability plot
         loglog          = FALSE,                 # log(-log Survival) plot
         logt            = FALSE,                 # log time
         time.inc        = 10,                   # time increment
         dots            = FALSE,                 # dot grid
         n.risk          = TRUE,                  # show number at risk
         ## srt.n.risk   = 0, sep.n.risk = 0.056, adj.n.risk = 1,
         ## y.n.risk     = 0, cex.n.risk = 0.6
         )

回帰分析

線形回帰モデル(Ordinary least square linear regression)

通常の線形回帰モデルはlm()関数を用います。

## 1982年の体重をいくつかの変数で回帰分析します。
lm1 <- lm(wt82 ~ wt71 + ht + age + sex + race, data = nhefs)

## 結果オブジェクトの名前だけ評価すると係数のみが出ます
lm1

## 
## Call:
## lm(formula = wt82 ~ wt71 + ht + age + sex + race, data = nhefs)
## 
## Coefficients:
## (Intercept)         wt71           ht          age          sex  
##     2.10715      0.88615      0.09027     -0.14217     -0.78569  
##        race  
##     0.30525

## summaryを使って詳しい結果を表示します。
## Wald testによるp値が表示されます。
summary(lm1)

## 
## Call:
## lm(formula = wt82 ~ wt71 + ht + age + sex + race, data = nhefs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.142  -4.218  -0.372   3.872  45.869 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.10715    5.26830   0.400  0.68923    
## wt71         0.88615    0.01385  63.969  < 2e-16 ***
## ht           0.09027    0.03038   2.971  0.00301 ** 
## age         -0.14217    0.01553  -9.158  < 2e-16 ***
## sex         -0.78569    0.52483  -1.497  0.13457    
## race         0.30525    0.53784   0.568  0.57041    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.439 on 1673 degrees of freedom
##   (67 observations deleted due to missingness)
## Multiple R-squared:  0.7871, Adjusted R-squared:  0.7865 
## F-statistic:  1237 on 5 and 1673 DF,  p-value: < 2.2e-16

## 信頼区間はconfintを使用します。
confint(lm1)

##                   2.5 %     97.5 %
## (Intercept) -8.22599760 12.4403034
## wt71         0.85898151  0.9133230
## ht           0.03067666  0.1498550
## age         -0.17262488 -0.1117225
## sex         -1.81508913  0.2437029
## race        -0.74965019  1.3601526

## tableoneがわかりやすいかもしれません。
library(tableone)
ShowRegTable(lm1, exp = FALSE)

##             beta [confint]       p     
## (Intercept)  2.11 [-8.23, 12.44]  0.689
## wt71         0.89 [0.86, 0.91]   <0.001
## ht           0.09 [0.03, 0.15]    0.003
## age         -0.14 [-0.17, -0.11] <0.001
## sex         -0.79 [-1.82, 0.24]   0.135
## race         0.31 [-0.75, 1.36]   0.570

## 結果オブジェクトはリストになっています。何が含まれているかチェックします。
names(lm1)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "na.action"     "xlevels"       "call"          "terms"        
## [13] "model"

## 各患者への予測値を表示します。
head(lm1$fitted.values)

##        1        2        3        4        5        6 
## 82.20576 63.33012 59.21709 60.76188 90.01230 97.88231

## summaryの結果も同様にリストです。
names(summary(lm1))

##  [1] "call"          "terms"         "residuals"     "coefficients" 
##  [5] "aliased"       "sigma"         "df"            "r.squared"    
##  [9] "adj.r.squared" "fstatistic"    "cov.unscaled"  "na.action"

## 要素を取り出せます。
summary(lm1)$coefficients

##                Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)  2.10715291 5.26829828  0.3999684 6.892308e-01
## wt71         0.88615225 0.01385285 63.9689613 0.000000e+00
## ht           0.09026583 0.03038120  2.9711080 3.009589e-03
## age         -0.14217367 0.01552538 -9.1575026 1.504687e-19
## sex         -0.78569311 0.52483173 -1.4970381 1.345720e-01
## race         0.30525121 0.53783551  0.5675550 5.704133e-01

## 任意の患者に対して予測を行うにはpredictを使用します。
predict(lm1, newdata = data.frame(wt71 = 100, ht = 170, age = 35, sex = 1, race = 1))

##       1 
## 100.611

## 交互作用は : を使います。
lm2 <- lm(wt82 ~ wt71 + ht + age + sex + race + sex:race, data = nhefs)
summary(lm2)

## 
## Call:
## lm(formula = wt82 ~ wt71 + ht + age + sex + race + sex:race, 
##     data = nhefs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.265  -4.164  -0.415   3.869  45.795 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.94980    5.26757   0.370   0.7113    
## wt71         0.88408    0.01392  63.514   <2e-16 ***
## ht           0.09243    0.03041   3.040   0.0024 ** 
## age         -0.14121    0.01553  -9.090   <2e-16 ***
## sex         -0.99328    0.54337  -1.828   0.0677 .  
## race        -0.55141    0.79340  -0.695   0.4872    
## sex:race     1.58445    1.07913   1.468   0.1422    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.436 on 1672 degrees of freedom
##   (67 observations deleted due to missingness)
## Multiple R-squared:  0.7874, Adjusted R-squared:  0.7866 
## F-statistic:  1032 on 6 and 1672 DF,  p-value: < 2.2e-16

## このように簡略化できます.
lm3 <- lm(wt82 ~ wt71 + ht + age + sex*race, data = nhefs)
summary(lm3)

## 
## Call:
## lm(formula = wt82 ~ wt71 + ht + age + sex * race, data = nhefs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.265  -4.164  -0.415   3.869  45.795 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.94980    5.26757   0.370   0.7113    
## wt71         0.88408    0.01392  63.514   <2e-16 ***
## ht           0.09243    0.03041   3.040   0.0024 ** 
## age         -0.14121    0.01553  -9.090   <2e-16 ***
## sex         -0.99328    0.54337  -1.828   0.0677 .  
## race        -0.55141    0.79340  -0.695   0.4872    
## sex:race     1.58445    1.07913   1.468   0.1422    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.436 on 1672 degrees of freedom
##   (67 observations deleted due to missingness)
## Multiple R-squared:  0.7874, Adjusted R-squared:  0.7866 
## F-statistic:  1032 on 6 and 1672 DF,  p-value: < 2.2e-16

## ネストしているモデルの比較にはanova()を使用します。
## この場合は交互作用項を検定しています。
anova(lm1, lm3)

## Analysis of Variance Table
## 
## Model 1: wt82 ~ wt71 + ht + age + sex + race
## Model 2: wt82 ~ wt71 + ht + age + sex * race
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1   1673 92578                           
## 2   1672 92459  1    119.21 2.1558 0.1422

## カテゴリカル変数はあらかじめfactorにしておくか、その場で変換します。
## MARITAL STATUS IN 1971  1: Under 17, 2: Married, 3: Widowed, 4: Never married, 5: Divorced, 6: Separated, 8: Unknown
lm4 <- lm(wt82 ~ wt71 + ht + age + sex + race + factor(marital), data = nhefs)
summary(lm4)

## 
## Call:
## lm(formula = wt82 ~ wt71 + ht + age + sex + race + factor(marital), 
##     data = nhefs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.820  -4.112  -0.343   3.935  45.913 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       2.10190    5.27285   0.399   0.6902    
## wt71              0.88644    0.01387  63.913   <2e-16 ***
## ht                0.09065    0.03040   2.982   0.0029 ** 
## age              -0.14530    0.01622  -8.959   <2e-16 ***
## sex              -0.87082    0.52993  -1.643   0.1005    
## race              0.39156    0.55203   0.709   0.4782    
## factor(marital)3  0.72344    0.84010   0.861   0.3893    
## factor(marital)4  0.08283    0.77109   0.107   0.9145    
## factor(marital)5  1.33100    0.77386   1.720   0.0856 .  
## factor(marital)6 -1.07020    1.01703  -1.052   0.2928    
## factor(marital)8 -0.94285    7.45122  -0.127   0.8993    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.439 on 1668 degrees of freedom
##   (67 observations deleted due to missingness)
## Multiple R-squared:  0.7877, Adjusted R-squared:  0.7865 
## F-statistic:   619 on 10 and 1668 DF,  p-value: < 2.2e-16

## カテゴリカル変数全体を検定したい場合はdrop1を使用できます。
## degrees of freedomが5になっているのがわかります。
drop1(lm4, test = "F")

## Single term deletions
## 
## Model:
## wt82 ~ wt71 + ht + age + sex + race + factor(marital)
##                 Df Sum of Sq    RSS    AIC   F value    Pr(>F)    
## <none>                        92301 6749.5                        
## wt71             1    226042 318343 8826.2 4084.8504 < 2.2e-16 ***
## ht               1       492  92793 6756.4    8.8934  0.002904 ** 
## age              1      4442  96743 6826.4   80.2705 < 2.2e-16 ***
## sex              1       149  92451 6750.2    2.7004  0.100515    
## race             1        28  92329 6748.0    0.5031  0.478229    
## factor(marital)  5       277  92578 6744.5    1.0005  0.415952    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## anova()でも同様です。
anova(lm1, lm4)

## Analysis of Variance Table
## 
## Model 1: wt82 ~ wt71 + ht + age + sex + race
## Model 2: wt82 ~ wt71 + ht + age + sex + race + factor(marital)
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1   1673 92578                           
## 2   1668 92301  5    276.82 1.0005  0.416

## plot()でモデル評価を表示できます。4つまとめて表示します。
layout(matrix(1:4, ncol = 2))
plot(lm1)

より詳しいモデル評価は、こちらの使用を参考にしてください http://rpubs.com/kaz_yos/residuals-bio213

一般化線形回帰モデル(Generalized linear models)

一般化線形モデルはglm()関数で行います。ロジスティック回帰分析やポアソン回帰分析は一般化線形モデルの枠組みで行います。

ロジスティック回帰分析

## 死亡の有無の二値のアウトカムをロジスティック回帰分析します。
glm1 <- glm(death ~ wt71 + ht + age + sex + race, data = nhefs,
            family = binomial(link = "logit"))
summary(glm1)

## 
## Call:
## glm(formula = death ~ wt71 + ht + age + sex + race, family = binomial(link = "logit"), 
##     data = nhefs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6940  -0.5883  -0.3127  -0.1699   2.9626  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.381923   2.080002  -2.107 0.035144 *  
## wt71         0.008342   0.005218   1.599 0.109884    
## ht          -0.018211   0.011832  -1.539 0.123753    
## age          0.118659   0.007379  16.081  < 2e-16 ***
## sex         -0.788801   0.209252  -3.770 0.000163 ***
## race         0.395281   0.203184   1.945 0.051723 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1707.3  on 1745  degrees of freedom
## Residual deviance: 1272.9  on 1740  degrees of freedom
## AIC: 1284.9
## 
## Number of Fisher Scoring iterations: 5

## 信頼区間はconfintです。
## これは通常の beta +/- 1.96*SE(beta)のWald信頼区間ではなく
## likelihood ratio testに対応するprofile likelihood信頼区間です
## このためsummary()のp値(Wald p)と食い違う場合があります
confint(glm1)

## Waiting for profiling to be done...

##                    2.5 %       97.5 %
## (Intercept) -8.473918006 -0.314061755
## wt71        -0.002019199  0.018451767
## ht          -0.041498208  0.004916647
## age          0.104557032  0.133509189
## sex         -1.201246484 -0.380320154
## race        -0.008416247  0.789329414

## tableoneもprofile likelihoodを用います。
ShowRegTable(glm1)

##             exp(beta) [confint] p     
## (Intercept) 0.01 [0.00, 0.73]    0.035
## wt71        1.01 [1.00, 1.02]    0.110
## ht          0.98 [0.96, 1.00]    0.124
## age         1.13 [1.11, 1.14]   <0.001
## sex         0.45 [0.30, 0.68]   <0.001
## race        1.48 [0.99, 2.20]    0.052

## Wald信頼区間はconfint.default()で得られます。
confint.default(glm1)

##                    2.5 %       97.5 %
## (Intercept) -8.458653109 -0.305193160
## wt71        -0.001885006  0.018569339
## ht          -0.041401131  0.004978226
## age          0.104196500  0.133121617
## sex         -1.198927806 -0.378674541
## race        -0.002952405  0.793514655

## likelihood ratio testのpはprofile likelihood信頼区間と一致します。
drop1(glm1, test = "LRT")

## Single term deletions
## 
## Model:
## death ~ wt71 + ht + age + sex + race
##        Df Deviance    AIC    LRT  Pr(>Chi)    
## <none>      1273.0 1285.0                     
## wt71    1   1275.5 1285.5   2.50 0.1136191    
## ht      1   1275.3 1285.3   2.38 0.1229251    
## age     1   1652.8 1662.8 379.90 < 2.2e-16 ***
## sex     1   1287.3 1297.3  14.39 0.0001483 ***
## race    1   1276.6 1286.6   3.69 0.0548903 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ポアソン回帰分析

## Introductory Statistics with Rパッケージ
library(ISwR)

## 
## Attaching package: 'ISwR'
## 
## The following object is masked from 'package:survival':
## 
##     lung

## デンマークの肺癌のデータを読み込みます。
## popがperson-timeと考えられますので、case count/person-timeのrate dataです。
data(eba1977)
eba1977

##          city   age  pop cases
## 1  Fredericia 40-54 3059    11
## 2     Horsens 40-54 2879    13
## 3     Kolding 40-54 3142     4
## 4       Vejle 40-54 2520     5
## 5  Fredericia 55-59  800    11
## 6     Horsens 55-59 1083     6
## 7     Kolding 55-59 1050     8
## 8       Vejle 55-59  878     7
## 9  Fredericia 60-64  710    11
## 10    Horsens 60-64  923    15
## 11    Kolding 60-64  895     7
## 12      Vejle 60-64  839    10
## 13 Fredericia 65-69  581    10
## 14    Horsens 65-69  834    10
## 15    Kolding 65-69  702    11
## 16      Vejle 65-69  631    14
## 17 Fredericia 70-74  509    11
## 18    Horsens 70-74  634    12
## 19    Kolding 70-74  535     9
## 20      Vejle 70-74  539     8
## 21 Fredericia   75+  605    10
## 22    Horsens   75+  782     2
## 23    Kolding   75+  659    12
## 24      Vejle   75+  619     7

## person-timeはlog変換してoffset項として取り込みます。
glm2 <- glm(formula = cases ~ city + age + offset(log(pop)),
            family  = poisson(link = "log"),
            data    = eba1977)
summary(glm2)

## 
## Call:
## glm(formula = cases ~ city + age + offset(log(pop)), family = poisson(link = "log"), 
##     data = eba1977)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.63573  -0.67296  -0.03436   0.37258   1.85267  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.6321     0.2003 -28.125  < 2e-16 ***
## cityHorsens  -0.3301     0.1815  -1.818   0.0690 .  
## cityKolding  -0.3715     0.1878  -1.978   0.0479 *  
## cityVejle    -0.2723     0.1879  -1.450   0.1472    
## age55-59      1.1010     0.2483   4.434 9.23e-06 ***
## age60-64      1.5186     0.2316   6.556 5.53e-11 ***
## age65-69      1.7677     0.2294   7.704 1.31e-14 ***
## age70-74      1.8569     0.2353   7.891 3.00e-15 ***
## age75+        1.4197     0.2503   5.672 1.41e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 129.908  on 23  degrees of freedom
## Residual deviance:  23.447  on 15  degrees of freedom
## AIC: 137.84
## 
## Number of Fisher Scoring iterations: 5

## 40-54歳のreferenceと比べて肺癌のlog rate ratioは年齢が上ると上るようです。
ShowRegTable(glm2)

##             exp(beta) [confint] p     
## (Intercept) 0.00 [0.00, 0.01]   <0.001
## cityHorsens 0.72 [0.50, 1.03]    0.069
## cityKolding 0.69 [0.48, 1.00]    0.048
## cityVejle   0.76 [0.53, 1.10]    0.147
## age55-59    3.01 [1.84, 4.90]   <0.001
## age60-64    4.57 [2.91, 7.24]   <0.001
## age65-69    5.86 [3.75, 9.25]   <0.001
## age70-74    6.40 [4.04, 10.21]  <0.001
## age75+      4.14 [2.52, 6.76]   <0.001

Risk regression

binomial(link = “log”)で理論的にはrisk ratioが出るrisk regressionが出来るはずですが、多くの場合アルゴリズムが収束しないです。このためポアソン回帰分析をしてrobust variance estimatorを組合せるmodified Poisson regression(Zho et al)が行われるようです。

## 二値アウトカムにたいしてポアソン回帰分析
glm3 <- glm(death ~ wt71 + ht + age + sex + race + factor(marital), data = nhefs,
            family = poisson(link = "log"))
summary(glm3)

## 
## Call:
## glm(formula = death ~ wt71 + ht + age + sex + race + factor(marital), 
##     family = poisson(link = "log"), data = nhefs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5974  -0.5432  -0.3543  -0.2358   2.3618  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -4.292221   1.605534  -2.673  0.00751 ** 
## wt71               0.006634   0.004086   1.623  0.10450    
## ht                -0.009480   0.009071  -1.045  0.29598    
## age                0.078774   0.005285  14.906  < 2e-16 ***
## sex               -0.473735   0.164039  -2.888  0.00388 ** 
## race               0.179196   0.156787   1.143  0.25307    
## factor(marital)3   0.039501   0.183079   0.216  0.82917    
## factor(marital)4   0.004181   0.285905   0.015  0.98833    
## factor(marital)5   0.363274   0.223957   1.622  0.10479    
## factor(marital)6   0.468670   0.283181   1.655  0.09792 .  
## factor(marital)8  -9.367165 469.323654  -0.020  0.98408    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1106.14  on 1745  degrees of freedom
## Residual deviance:  783.59  on 1735  degrees of freedom
## AIC: 1475.6
## 
## Number of Fisher Scoring iterations: 11

## robust sandwich covariance matrix estimatorのパッケージ
library(sandwich)
## 任意のvariance estimatorを用いて検定するパッケージ
library(lmtest)

## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## 
## Attaching package: 'lmtest'
## 
## The following object is masked from 'package:rms':
## 
##     lrtest

## 実際に検定する。標準誤差が小さくなりました。
## Estimateはlog risk ratioとして解釈されます。
lmtest1 <- coeftest(glm3, vcov = sandwich)
lmtest1

## 
## z test of coefficients:
## 
##                    Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)      -4.2922208  1.2277391 -3.4960 0.0004722 ***
## wt71              0.0066337  0.0038200  1.7366 0.0824614 .  
## ht               -0.0094801  0.0069847 -1.3573 0.1746980    
## age               0.0787744  0.0040701 19.3543 < 2.2e-16 ***
## sex              -0.4737350  0.1290475 -3.6710 0.0002416 ***
## race              0.1791961  0.1201800  1.4911 0.1359446    
## factor(marital)3  0.0395015  0.1311076  0.3013 0.7631931    
## factor(marital)4  0.0041811  0.2544032  0.0164 0.9868873    
## factor(marital)5  0.3632743  0.1943636  1.8690 0.0616166 .  
## factor(marital)6  0.4686700  0.1989710  2.3555 0.0184993 *  
## factor(marital)8 -9.3671652  1.0116367 -9.2594 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## exp変換と信頼区間は自分で計算する必要があります。
cbind(exp(lmtest1[,1]),
      ## lower
      exp(lmtest1[,1] - 1.96*lmtest1[,2]),
      ## upper
      exp(lmtest1[,1] + 1.96*lmtest1[,2]))

##                           [,1]          [,2]         [,3]
## (Intercept)      0.01367452363 0.00123264946 0.1516997350
## wt71             1.00665575048 0.99914690157 1.0142210303
## ht               0.99056467508 0.97709616029 1.0042188429
## age              1.08196025613 1.07336333133 1.0906260366
## sex              0.62267225736 0.48351739236 0.8018754779
## race             1.19625533918 0.94520195530 1.5139905588
## factor(marital)3 1.04029200904 0.80455217904 1.3451053795
## factor(marital)4 1.00418989703 0.60990680469 1.6533630082
## factor(marital)5 1.43803026197 0.98247690476 2.1048138886
## factor(marital)6 1.59786758395 1.08186550716 2.3599798672
## factor(marital)8 0.00008548538 0.00001176981 0.0006208896

生存分析

Cox回帰分析はsurvivalパッケージのcoxph()を用います。

library(survival)

## 大腸癌の生存データを用います。
data(colon)

## Surv(time, event_status) ~ という書きかたをします
coxph1 <- coxph(Surv(time, status) ~ rx + sex + age + obstruct + perfor, data = colon)
summary(coxph1)

## Call:
## coxph(formula = Surv(time, status) ~ rx + sex + age + obstruct + 
##     perfor, data = colon)
## 
##   n= 1858, number of events= 920 
## 
##                coef exp(coef)  se(coef)      z    Pr(>|z|)    
## rxLev     -0.019291  0.980894  0.076887 -0.251     0.80190    
## rxLev+5FU -0.437009  0.645966  0.083987 -5.203 0.000000196 ***
## sex       -0.043975  0.956978  0.066199 -0.664     0.50651    
## age       -0.001236  0.998764  0.002815 -0.439     0.66048    
## obstruct   0.211936  1.236068  0.081680  2.595     0.00947 ** 
## perfor     0.198608  1.219703  0.181990  1.091     0.27514    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##           exp(coef) exp(-coef) lower .95 upper .95
## rxLev        0.9809     1.0195    0.8437    1.1404
## rxLev+5FU    0.6460     1.5481    0.5479    0.7616
## sex          0.9570     1.0450    0.8405    1.0896
## age          0.9988     1.0012    0.9933    1.0043
## obstruct     1.2361     0.8090    1.0532    1.4507
## perfor       1.2197     0.8199    0.8538    1.7425
## 
## Concordance= 0.563  (se = 0.01 )
## Rsquare= 0.024   (max possible= 0.999 )
## Likelihood ratio test= 44.77  on 6 df,   p=0.00000005196
## Wald test            = 43.17  on 6 df,   p=0.0000001078
## Score (logrank) test = 43.75  on 6 df,   p=0.00000008292

## 比例ハザード性の検定が出きます
cox.zph(coxph1)

##                rho   chisq       p
## rxLev     -0.02779  0.7111 0.39908
## rxLev+5FU -0.01935  0.3450 0.55697
## sex        0.06289  3.6588 0.05577
## age       -0.00133  0.0018 0.96612
## obstruct  -0.08963  7.6356 0.00572
## perfor     0.02179  0.4543 0.50032
## GLOBAL          NA 12.6781 0.04844

## このオブジェクトをplotするとlog hazard ratioの時間依存性を確認できます。
## plot(cox.zph(coxph1))

## 連続変数をより柔軟にするにはnatural splineが使えます。
library(splines)
coxph2 <- coxph(Surv(time, status) ~ rx + sex + ns(age, df = 3) + obstruct + perfor, data = colon)
summary(coxph2)

## Call:
## coxph(formula = Surv(time, status) ~ rx + sex + ns(age, df = 3) + 
##     obstruct + perfor, data = colon)
## 
##   n= 1858, number of events= 920 
## 
##                      coef exp(coef) se(coef)      z    Pr(>|z|)    
## rxLev            -0.01122   0.98885  0.07699 -0.146     0.88416    
## rxLev+5FU        -0.44218   0.64263  0.08402 -5.263 0.000000142 ***
## sex              -0.03784   0.96286  0.06635 -0.570     0.56842    
## ns(age, df = 3)1 -0.26691   0.76574  0.15832 -1.686     0.09181 .  
## ns(age, df = 3)2 -1.70749   0.18132  0.62423 -2.735     0.00623 ** 
## ns(age, df = 3)3 -0.23918   0.78727  0.24120 -0.992     0.32139    
## obstruct          0.20464   1.22708  0.08174  2.504     0.01230 *  
## perfor            0.19848   1.21954  0.18192  1.091     0.27528    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## rxLev               0.9888     1.0113   0.85034    1.1499
## rxLev+5FU           0.6426     1.5561   0.54506    0.7577
## sex                 0.9629     1.0386   0.84545    1.0966
## ns(age, df = 3)1    0.7657     1.3059   0.56146    1.0443
## ns(age, df = 3)2    0.1813     5.5151   0.05335    0.6163
## ns(age, df = 3)3    0.7873     1.2702   0.49070    1.2631
## obstruct            1.2271     0.8149   1.04543    1.4403
## perfor              1.2195     0.8200   0.85378    1.7420
## 
## Concordance= 0.571  (se = 0.01 )
## Rsquare= 0.028   (max possible= 0.999 )
## Likelihood ratio test= 52.79  on 8 df,   p=0.00000001186
## Wald test            = 52.22  on 8 df,   p=0.00000001527
## Score (logrank) test = 52.77  on 8 df,   p=0.00000001193

## splineをグラフにするにはexpand.gridで年齢を変化させたデータセットを作成
newdat <- expand.grid(age = seq(from = 18, to = 85, by = 1),
                      rx = "Lev", sex = 1, obstruct = 0, perfor = 0)
## predictでアウトカムを予測riskとありますがhazard ratioのことのようです。
## くわしくは ?predict.coxph
newdat$risk <- predict(coxph2, newdat, type = "risk")
## plotします。若年発症の大腸癌は大変予後が悪いということでしょうか。
plot(risk ~ age, data =  newdat)

## 二次項を入れるだけならpolyが使用できます。
## collinearityを起こさないように変換した変数が投入されます。
coxph3 <- coxph(Surv(time, status) ~ rx + sex + poly(age, 2) + obstruct + perfor, data = colon)
summary(coxph3)

## Call:
## coxph(formula = Surv(time, status) ~ rx + sex + poly(age, 2) + 
##     obstruct + perfor, data = colon)
## 
##   n= 1858, number of events= 920 
## 
##                   coef exp(coef) se(coef)      z    Pr(>|z|)    
## rxLev         -0.01511   0.98500  0.07693 -0.196      0.8443    
## rxLev+5FU     -0.43936   0.64445  0.08400 -5.230 0.000000169 ***
## sex           -0.03627   0.96438  0.06632 -0.547      0.5844    
## poly(age, 2)1 -0.57079   0.56508  1.38558 -0.412      0.6804    
## poly(age, 2)2  3.20288  24.60329  1.35459  2.364      0.0181 *  
## obstruct       0.20596   1.22870  0.08174  2.520      0.0117 *  
## perfor         0.20161   1.22337  0.18192  1.108      0.2678    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## rxLev            0.9850    1.01522   0.84713    1.1453
## rxLev+5FU        0.6445    1.55171   0.54662    0.7598
## sex              0.9644    1.03694   0.84683    1.0982
## poly(age, 2)1    0.5651    1.76966   0.03739    8.5411
## poly(age, 2)2   24.6033    0.04064   1.72966  349.9651
## obstruct         1.2287    0.81387   1.04682    1.4422
## perfor           1.2234    0.81741   0.85646    1.7475
## 
## Concordance= 0.57  (se = 0.01 )
## Rsquare= 0.027   (max possible= 0.999 )
## Likelihood ratio test= 50.11  on 7 df,   p=0.00000001375
## Wald test            = 49.02  on 7 df,   p=0.00000002248
## Score (logrank) test = 49.59  on 7 df,   p=0.00000001737

newdat$risk2 <- predict(coxph3, newdat, type = "risk")
plot(risk2 ~ age, data =  newdat)

線形混合効果モデル(Linear mixed effects model)

連続変数のアウトカムの場合はlme4パッケージが使用可能です。GLMMの場合はMCMCglmmパッケージや外部ソフトウェアのStanなどMCMC samplerを用いた方が良いと思われます。

## Class roomデータ: http://www-personal.umich.edu/~bwest/chapter4.html
## school/class/studentの三層構造の数学の成績向上のデータ
classroom <- read.csv("http://www-personal.umich.edu/~bwest/classroom.csv")

## 必要なパッケージの読み込み
library(lme4)

## Loading required package: Matrix
## Loading required package: Rcpp

library(lmerTest)

## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step

## Random intercept model
## Random effectsを見るとアウトカムのばらつきの分布が見られます。
## 数学の点数の向上のクラスによるばらつき(classid)と学校によるばらつき(schoolid)は同じぐらい
## 生徒個々人のばらつき(Residualがlevel 1のばらつき)の方がずっと大きいようです。
lmer1 <- lmer(mathgain ~ mathkind + sex + minority + ses + housepov + (1 | schoolid/classid),
              data = classroom)
summary(lmer1)

## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [merModLmerTest]
## Formula: mathgain ~ mathkind + sex + minority + ses + housepov + (1 |  
##     schoolid/classid)
##    Data: classroom
## 
## REML criterion at convergence: 11378.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.7926 -0.6052 -0.0327  0.5524  4.2727 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  classid:schoolid (Intercept)  81.56    9.031  
##  schoolid         (Intercept)  77.76    8.818  
##  Residual                     734.42   27.100  
## Number of obs: 1190, groups:  classid:schoolid, 312; schoolid, 107
## 
## Fixed effects:
##               Estimate Std. Error         df t value  Pr(>|t|)    
## (Intercept)  285.05797   11.02077 1138.70000  25.866   < 2e-16 ***
## mathkind      -0.47086    0.02228 1173.10000 -21.133   < 2e-16 ***
## sex           -1.23459    1.65743 1126.20000  -0.745   0.45650    
## minority      -7.75587    2.38499  643.70000  -3.252   0.00121 ** 
## ses            5.23971    1.24497 1171.20000   4.209 0.0000276 ***
## housepov     -11.43920    9.93736  119.50000  -1.151   0.25198    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) mthknd sex    minrty ses   
## mathkind -0.969                            
## sex      -0.042 -0.032                     
## minority -0.265  0.153 -0.015              
## ses       0.123 -0.165  0.019  0.144       
## housepov -0.173  0.035 -0.009 -0.184  0.078

## random intercept/slope model
## Minorityであることの影響は平均として-7.8点ぐらいのようです。このminorityであることの影響の
## 傾きのばらつきはschool level > class levelのようです。学校によってminorityの割り合いが違う
## とかそういった生徒への教育方針が違ったりするのかも?
lmer2 <- lmer(mathgain ~ mathkind + sex + minority + ses + housepov + (1 + minority | schoolid/classid),
              data = classroom)
summary(lmer2)

## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [merModLmerTest]
## Formula: mathgain ~ mathkind + sex + minority + ses + housepov + (1 +  
##     minority | schoolid/classid)
##    Data: classroom
## 
## REML criterion at convergence: 11375
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.7763 -0.6090 -0.0108  0.5526  4.2747 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev. Corr 
##  classid:schoolid (Intercept)  62.58    7.911        
##                   minority     28.82    5.369   0.07 
##  schoolid         (Intercept) 119.12   10.914        
##                   minority    111.45   10.557   -0.68
##  Residual                     721.48   26.860        
## Number of obs: 1190, groups:  classid:schoolid, 312; schoolid, 107
## 
## Fixed effects:
##               Estimate Std. Error         df t value  Pr(>|t|)    
## (Intercept)  285.94776   11.05968 1086.00000  25.855   < 2e-16 ***
## mathkind      -0.47296    0.02228 1158.90000 -21.227   < 2e-16 ***
## sex           -1.20077    1.65464 1119.00000  -0.726   0.46817    
## minority      -7.88572    2.72755   62.70000  -2.891   0.00527 ** 
## ses            5.04625    1.24562 1160.00000   4.051 0.0000543 ***
## housepov     -10.87370    9.95607  112.70000  -1.092   0.27709    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) mthknd sex    minrty ses   
## mathkind -0.966                            
## sex      -0.037 -0.035                     
## minority -0.277  0.139 -0.019              
## ses       0.123 -0.163  0.020  0.120       
## housepov -0.163  0.033 -0.012 -0.193  0.079

一般化推定方程式(Generalized estimating equation)

geepackパッケージが使用可能です。

library(geepack)

## Handbook of Statistical Analyses Using R 2nd ed. この本はインターネットで見られます。
library(HSAUR2)

## Loading required package: tools

## 詳細不明ですが呼吸器疾患のplacebo controlled trial、アウトカムは二値です。
data(respiratory, package = "HSAUR2")

## 数値に変換しておきます
respiratory$status <- as.numeric(respiratory$status == "good")
respiratory$month <- as.numeric(as.character(respiratory$month))

## 時間の影響がgood statusのlog oddsに直線的に影響を与えるとします。
## 治療と時間の変数の交互作用は有意で両群の経過が異なることを示しているようです。
geeglm1 <- geeglm(formula   = status ~ treatment*month + gender,
                  family    = binomial(link = "logit"),
                  id        = subject,
                  data      = respiratory,
                  corstr    = "exchangeable")
summary(geeglm1)

## 
## Call:
## geeglm(formula = status ~ treatment * month + gender, family = binomial(link = "logit"), 
##     data = respiratory, id = subject, corstr = "exchangeable")
## 
##  Coefficients:
##                          Estimate  Std.err  Wald Pr(>|W|)  
## (Intercept)              -0.12631  0.24789 0.260   0.6104  
## treatmenttreatment        0.36685  0.33498 1.199   0.2735  
## month                    -0.02846  0.06591 0.187   0.6658  
## gendermale               -0.11855  0.35977 0.109   0.7418  
## treatmenttreatment:month  0.20640  0.10060 4.210   0.0402 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Estimated Scale Parameters:
##             Estimate Std.err
## (Intercept)    1.004 0.04357
## 
## Correlation: Structure = exchangeable  Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha   0.4671 0.05675
## Number of clusters:   111   Maximum cluster size: 5

グラフ作成

利用可能なグラフシステムの種類

Rにはグラフ作成のためのシステムがいくつかあります。おおまかにいって、traditional (base)と呼ばれるシステムとgridと呼ばれるシステムがあります。gridの方は最近は、実質的にはggplot2というパッケージが主に使われます。

baseグラフィックによる作図

通常のplot関数はbaseの方になります。formulaを使ったグラフ指定ができます。plot関数はオブジェクトの種類に応じてさまざまなものをplotできます。統計モデルに使用するとモデル評価を表示しますがこれらはモデリングのところで扱います。先程すでに出たカプランマイヤー曲線の描画もbaseを使用しています。

## 連続変数 ~ カテゴリカル変数は自動的にboxplotになります。
plot(time ~ x, data = aml)

## 連続変数同士は散布図になります
plot(wt71 ~ age, nhefs)

## 軸のラベルは下記のようにおこないます。
plot(wt71 ~ age, nhefs,
     xlab = "Age in 1971", ylab = "Weight in 1971", main = "Scatterplot", sub = "subtitle")

## 軸の範囲は下記のように指定します。
plot(wt71 ~ age, nhefs,
     xlim = c(0,100), ylim = c(0,200))

ggplot2パッケージによる作図

Grammer of Graphics (gg)というグラフをどのように言語として表現すべきかという考えかたがあります。これをRパッケージとして具現化したのがggplot2パッケージです。かなりくせがあるように感じますが、実に理路整然とした形でグラフを言語表現します。簡略記法もあるのですが、逆にggの概念を隠してしまうので、あえて完全記法で記載します。

## ggplot2パッケージを読み込みます。
library(ggplot2)

## デフォルトのデータセットとどの軸(scale)にどの変数を割り振るかを指定します。
## +サインで表現を繋いでいきます。
## レイヤーとして形態(geometry)が散布図(point)のものを加えます。
## baseと違ってグラフをオブジェクトとして保存できます。
plot1 <- ggplot(data = nhefs, mapping = aes(x = age, y = wt71)) +
    layer(geom = "point")

## 名前だけ評価すれば表示されます。
plot1

## themeという要素をつかってグラフの仕上げを調節します。
## theme_bwというプリセットで背景を白にできます。
plot1 + theme_bw()

## パネルを変数で割るにはfacetという要素をつかいます。
## 1変数でわるならfacet_wrap
plot1 + facet_wrap( ~ sex)

## 縦横に展開するならfacet_gridを使います
## ラベリングはlabsを使います。
plot1 + facet_grid(sex ~ qsmk) + labs(title = "sex by qsmk", x = "Age in 1971", y = "Weight in 1971")

## 男女で色分けしてみます。カテゴリカルで与える必要があるのでfactorを使用します。
plot2 <- ggplot(data = nhefs, mapping = aes(x = age, y = wt71, color = factor(sex))) +
    layer(geom = "point") + labs(color = "sex")
plot2

## パネルを分けた方が見易いですね。
plot2 + facet_wrap( ~ sex)

3Dグラフ

いくつか3Dグラフを描画する方法もあります。

library(scatterplot3d)

## with関数を使うと表記を簡略化できます
with(nhefs,
     scatterplot3d(x = age, y = ht, z = wt71))

## あるいは
scatterplot3d(x = nhefs$age, y = nhefs$ht, z = nhefs$wt71)

実用性は微妙ですがインタラクティブな3Dグラフを作成することもできます。

##
library(rgl)

## 男女で青赤に塗り分けています
with(nhefs,
     plot3d(x = age, y = ht, z = wt71,
            type = "s", size = 1,
            col = c("blue","red")[sex + 1]))

グラフの出力

グラフの出力はファイル形式に応じた関数を用いて行います。

## PDFの場合は複数ページ対応しています。
pdf(file = "plots.pdf", width = 7, height = 7, family = "sans")
plot(time ~ x, data = aml)
plot(wt71 ~ age, nhefs)
plot1
## ファイルを閉じて完了というシグナルを送る必要があります。
dev.off()


pdf(file = "plot2.pdf", width = 7, height = 7, family = "sans")
## 3Dグラフの見易い角度を探すためにforループで複数角度を変えながら作成してみます
for (i in seq(0, 360, 10)) {
   scatterplot3d(x = nhefs$age, y = nhefs$ht, z = nhefs$wt71, angle = i, main = i)
}
dev.off()


## pngの場合(jpegは関数名以外同じ)は一枚づつです。
png(file = "plot1.png", width = 700, height = 700, family = "sans")
plot1
dev.off()

png(file = "plot2.png", width = 700, height = 700, family = "sans")
plot2
dev.off()