最終更新日 ： 2007.9.10
連続型変数とカテゴリカル変数を混在させる
p95のBOX6.1の解析例
> dat <- excel.w(3) # Excelからデータの読み込み
Read 90 items
# データフレームの作成
> Leprosy <- data.frame(TREATMT=as.factor(dat[,1]), BACBEF=dat[,2], BACAFTER=dat[,3])
> Leprosy
TREATMT BACBEF BACAFTER
1 1 7.464 7.789
2 1 8.551 0.000
3 1 10.774 4.742
（途中省略）
> attach(Leprosy) # データフレームLeprosyにパスを通す
> box6.1 <- glm(BACAFTER ~ BACBEF + TREATMT)
> summary(aov(box6.1)) # 分散分析表
Df Sum Sq Mean Sq F value Pr(>F)
BACBEF 1 587.48 587.48 40.1800 1.033e-06 ***
TREATMT 2 83.35 41.67 2.8502 0.076 .
Residuals 26 380.15 14.62
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> library(car) # carパッケージの呼び出し
> Anova(box6.1, test="F") # TypeIIの平方和
Anova Table (Type II tests)
Response: BACAFTER
SS Df F Pr(>F)
BACBEF 515.01 1 35.2238 2.912e-06 ***
TREATMT 83.35 2 2.8502 0.076 .
Residuals 380.15 26
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> summary(box6.1) # 係数表
Call:
glm(formula = BACAFTER ~ BACBEF + TREATMT)
Deviance Residuals:
Min 1Q Median 3Q Max
-6.5177 -1.7218 -0.1963 2.2941 9.3728
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.6026 1.8815 -0.852 0.4021
BACBEF 0.8831 0.1488 5.935 2.91e-06 ***
TREATMT2 0.8641 1.7548 0.492 0.6265
TREATMT3 3.9060 1.7326 2.254 0.0328 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 14.62122)
Null deviance: 1050.98 on 29 degrees of freedom
Residual deviance: 380.15 on 26 degrees of freedom
AIC: 171.32
Number of Fisher Scoring iterations: 2
> detach(Leprosy) # パスの解除
p96の図6.1の作図例
> attach(Leprosy)
# plot()を用いて処置1のデータをプロットする
> plot(BACBEF[TREATMT==1], BACAFTER[TREATMT==1],
+ xlim=c(0,25), ylim=c(0,20),
+ xlab="BACBEF", ylab="BACAFTER",
+ pch=18)
# points()を用いて処置2のデータを図に追加
> points(BACBEF[TREATMT==2], BACAFTER[TREATMT==2], pch=17)
# 同じく処置3のデータをプロット
> points(BACBEF[TREATMT==3], BACAFTER[TREATMT==3], pch=8)
> x <- seq(0, 25, 0.1) # 直線の描写につかうx軸の値
> model1 <- -1.6026 + 0.8831*x + 0 # モデル1のy軸の値
> points(x, model1, type="l") # モデル1の直線を描く
> model2 <- -1.6026 + 0.8831*x + 0.8641
> points(x, model2, type="l", lty=2)
> model3 <- -1.6026 + 0.8831*x + 3.9060
> points(x, model3, type="l", lty=3)
> # 直線の凡例
> legend(17, 5, # 凡例を書く座標
+ c("model1", "model2", "model3"), # 凡例のラベル
+ lty=c(1, 2, 3), # 凡例のラインのスタイル
+ )
> # 点の凡例
> legend(0, 20,
+ c("treatmt1", "treatmt2", "treatmt3"),
+ pch=c(18, 17, 8) # 点のスタイル
+ )
> detach(Reprosy)
p98のBOX 6.2の解析例および図の描写
# データフレームの作成
> fats <- data.frame(WEIGHT=dat[,1], FAT=dat[,2], SEX=as.factor(dat[,3]))
> fats
WEIGHT FAT SEX
1 89 28 2
2 88 27 2
3 66 24 2
4 59 23 2
（途中省略）
> attach(fats) # データフレームfatsにパスを通す
> box6.2 <- glm(FAT ~ WEIGHT + SEX) # 一般化線形モデルとして解析
> summary(aov(box6.2)) # 分散分析表の出力
Df Sum Sq Mean Sq F value Pr(>F)
WEIGHT 1 1.328 1.328 0.5184 0.4819
SEX 1 176.098 176.098 68.7293 3.482e-07 ***
Residuals 16 40.995 2.562
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> library(car) # carパッケージの呼び出し
> Anova(box6.2, test="F") # type2の平方和
Anova Table (Type II tests)
Response: FAT
SS Df F Pr(>F)
WEIGHT 87.105 1 33.996 2.556e-05 ***
SEX 176.098 1 68.729 3.482e-07 ***
Residuals 40.995 16
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 男性のデータポイントをプロット
> plot(WEIGHT[SEX==1], FAT[SEX==1],
+ xlim=c(50,100), ylim=c(20,40),
+ xlab="体重", ylab="脂肪量", pch=18)
# 女性のデータポイントを追加
> points(WEIGHT[SEX==2], FAT[SEX==2], pch=15)
> summary(box6.2) #係数表の出力
Call:
glm(formula = FAT ~ WEIGHT + SEX)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.7740 -0.7760 -0.2532 0.8698 3.0546
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16.96228 2.41021 7.038 2.80e-06 ***
WEIGHT 0.21715 0.03724 5.831 2.56e-05 ***
SEX2 -7.90375 0.95337 -8.290 3.48e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 2.562192)
Null deviance: 218.421 on 18 degrees of freedom
Residual deviance: 40.995 on 16 degrees of freedom
AIC: 76.531
Number of Fisher Scoring iterations: 2
> x <- seq(50, 100, 0.1) # x軸のデータ
> men <- 16.96228 + 0.21715*x + 0 # 男性のy軸のデータ
> women <- 16.96228 + 0.21715*x - 7.90375 # 女性のy軸のデータ
> points(x, men, type="l") # 男性のデータポイントへの直線の当てはめ
> points(x, women, type="l") # 女性も同様
> detach(fats) # パスの解除
ここまで繰り返しRのglm()を使用して解析してきたので、そろそろ解析する変数が異なっても、それに対応してRでデータ解析をすることが可能であろう。そんなわけで、これ以降は新しい解析例や特に難しい（分かりにくい）解析例以外は各自でやってみてほしい。どのような解析であっても、glm()の使い方は基本的に同じであるから、今までの章についてのRでの解析例を読み返せば分かるはずである。したがって、練習問題も教科書の解答を見ながらやればRでも必ず解析できるので、練習問題の補足を書くことはしない。どうしても分からない部分があれば、このwebサイトに併設してある掲示板なりで聞いてもらえれば対応するつもりである。
