# ←この記号の後は次の改行までコメントになります。
# 線形回帰分析の実行
# 医学部学生の性別（Sex, 二項の名義変数）、入学成績（EntranceExam, 連続変数）、卒業
# 試験成績（GradExam, 連続変数）、医師国家試験の結果（NatnExam, 二項の名義変数）を
# 含む架空データを用います。
# ここから34行目まではサンプルデータセットの作成です。
rm(list=ls())
set.seed(1017)
Sex<-factor(runif(500)<0.6,labels=c("Female","Male"))
entseed<-
rnorm(n=500,mean=0.2,sd=1)+
I(Sex=="Male")*rnorm(n=500,mean=0.05,sd=0.5)+
I(Sex=="Female")*rnorm(n=500,mean=-0.05,sd=0.5)
EntranceExam.All<-round(100*exp(entseed)/(1+exp(entseed)),0)
Dataset.1<-
data.frame(
Sex=Sex[EntranceExam.All>=quantile(EntranceExam.All,probs=0.8)],
EntranceExam=EntranceExam.All[EntranceExam.All>=quantile(EntranceExam.All,probs=0.8)])
gradseed<-
rnorm(n=nrow(Dataset.1),mean=1,sd=0.4)+
(Dataset.1$EntranceExam-mean(Dataset.1$EntranceExam))/sd(Dataset.1$EntranceExam)*rnorm(n=nrow(Dataset.1),mean=0.1,sd=0.2)+
I(Dataset.1$Sex=="Female")*rnorm(n=nrow(Dataset.1),mean=0.2,sd=0.2)+
I(Dataset.1$Sex=="Male")*rnorm(n=nrow(Dataset.1),mean=-0.2,sd=0.2)
Dataset.1$GradExam<-round(100*exp(gradseed)/(1+exp(gradseed)),0)
natnseed<-
rnorm(n=nrow(Dataset.1),mean=1.5,sd=0.8)+
(Dataset.1$GradExam-mean(Dataset.1$GradExam))/sd(Dataset.1$GradExam)*rnorm(n=nrow(Dataset.1),mean=0.0,sd=0.4)
Dataset.1$NatnExam<-
factor(
round(100*exp(natnseed)/(1+exp(natnseed)),0)>=60,
labels=c("Failed","Passed")
)
# 入学試験成績の分布を見る。
# 最低点76点で、最低点付近が最頻値で点が上がるごとに概ね頻度が下がっていくのがわかります。
summary(Dataset.1$EntranceExam)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 76.00 79.00 84.00 84.16 88.00 95.00
hist(Dataset.1$EntranceExam)
# 卒業試験成績の分布を見る。
# こちらは71点付近に最頻値があり、ほぼ正規分布していそうです。
summary(Dataset.1$GradExam)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 47.00 66.00 71.00 71.07 77.00 95.00
hist(Dataset.1$GradExam)
# さて、あなたは医学部教員だったとして、入学試験成績と卒業試験成績には相関があるだろうか
# という疑問を覚え、解析してみようと思い立ちました。
# まずは、入学試験成績をX軸に卒業試験成績をY軸にして、xy軸のプロットをしてみましょう。
plot(GradExam~EntranceExam,data=Dataset.1)
# なんとなく右上がりにも見えますので、相関はあるのかもしれません。
# これを数値化してみます。
fit.lm.1<- # 結果を格納するオブジェクト名は任意に決めます
lm( # 線形回帰分析ではlm関数を使います。
GradExam~EntranceExam, # 卒業試験成績を応答変数に、入学試験成績を説明変数にします
data=Dataset.1 # データセットの指定
)
summary(fit.lm.1)
##
## Call:
## lm(formula = GradExam ~ EntranceExam, data = Dataset.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.8150 -5.0884 0.0361 5.8329 21.4617
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.5690 12.7621 3.492 0.000715 ***
## EntranceExam 0.3149 0.1513 2.081 0.040005 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.381 on 100 degrees of freedom
## Multiple R-squared: 0.0415, Adjusted R-squared: 0.03192
## F-statistic: 4.33 on 1 and 100 DF, p-value: 0.04
# 入学試験成績が1点良いごとに、卒業試験成績が0.3点良い有意な（P=0.040）相関関係が見られ
# ました。また、入学試験成績がゼロ点と仮定したときの卒業試験の点数は45点（y切片）でした。
# この線形回帰分析の結果を、abline関数を使って散布図にプロットすることもできます。
abline(fit.lm.1)
# 入学時のベースライン変数には性別（Sex）も含まれています。これも卒業試験の成績に影響す
# るかもしれません。仮説は女子のほうが成績が良いのでは？とします。
# まずは男女別にプロットしてみましょう。
plot(GradExam~EntranceExam, # 回帰式は同様ですが、
data=subset(Dataset.1,Sex=="Male"), # データは男性のみを指定します。
xlim=c(70,100), # グラフの描画範囲は、x軸が70−100を、
ylim=c(40,100), # y軸が40-100を指定します。
col=4) # 男女を区別するために男には青い色を着けます。
par(new=T) # plot関数は実行される度に新しいウィンドウに
# 描画されますが、このpar関数でこのオプション
# を指定することにより重ね描きをします。
plot(GradExam~EntranceExam, # 回帰式は同様ですが、
data=subset(Dataset.1,Sex=="Female"),# データは女性のみを指定します。
xlim=c(70,100), # グラフの描画範囲は、x軸が70−100を、
ylim=c(40,100), # y軸が40-100を指定します。
col=2, # 男女を区別するために女には赤い色を着けます
axes=F, # 男のグラフで軸とラベルはすでに描画されてい
xlab="",ylab="") # このプロットでは描きません。
# 仮説どおり、女子のほうが卒業試験成績が良さそうです。
# 再び回帰分析をしてみます。今回は多変量の重回帰分析ですが、文法はほぼ一緒です。
fit.lm.2<- # 結果を格納するオブジェクト名は任意に決めます
lm( # 線形回帰分析ではlm関数を使います。
GradExam~EntranceExam+Sex, # 卒業試験成績を応答変数に、性別と入学試験成績を説明変
# 数にします
data=Dataset.1 # データセットの指定
)
summary(fit.lm.2)
##
## Call:
## lm(formula = GradExam ~ EntranceExam + Sex, data = Dataset.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.5286 -6.0120 0.7869 6.2118 22.9880
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.2763 12.5035 3.941 0.000151 ***
## EntranceExam 0.2989 0.1469 2.034 0.044587 *
## SexMale -4.7630 1.7683 -2.694 0.008304 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.131 on 99 degrees of freedom
## Multiple R-squared: 0.1069, Adjusted R-squared: 0.08891
## F-statistic: 5.928 on 2 and 99 DF, p-value: 0.003701
# 性別もモデルに含めた後で、入学試験成績が1点良いごとに、卒業試験成績が0.3点良い
# 有意な（P=0.045）相関関係が見られました。性別は男性のほうが卒業試験成績が4.8点
# 低い（女性が4.8点高い）有意な（P=0.008）相関が見られました。
# 男女をそれぞれ解析するサブグループ解析も行えます。
# 男性
fit.lm.2.m<- # 結果を格納するオブジェクト名は任意に決めます
lm( # 線形回帰分析ではlm関数を使います。
GradExam~EntranceExam, # 卒業試験成績を応答変数に、入学試験成績を説明変数にします
data=subset(Dataset.1,Sex=="Male")
# 男性のデータセットの指定
)
summary(fit.lm.2.m)
##
## Call:
## lm(formula = GradExam ~ EntranceExam, data = subset(Dataset.1,
## Sex == "Male"))
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.0698 -6.2556 0.4657 6.7890 22.4657
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.0188 15.7052 2.484 0.0154 *
## EntranceExam 0.3643 0.1865 1.953 0.0548 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.787 on 70 degrees of freedom
## Multiple R-squared: 0.05168, Adjusted R-squared: 0.03813
## F-statistic: 3.814 on 1 and 70 DF, p-value: 0.05481
# 女性
fit.lm.2.f<- # 結果を格納するオブジェクト名は任意に決めます
lm( # 線形回帰分析ではlm関数を使います。
GradExam~EntranceExam, # 卒業試験成績を応答変数に、入学試験成績を説明変数にします
data=subset(Dataset.1,Sex=="Female")
# 女性のデータセットの指定
)
summary(fit.lm.2.f)
##
## Call:
## lm(formula = GradExam ~ EntranceExam, data = subset(Dataset.1,
## Sex == "Female"))
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.1042 -5.6882 -0.1528 5.0581 11.5153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63.8143 18.3208 3.483 0.00165 **
## EntranceExam 0.1269 0.2164 0.586 0.56242
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.285 on 28 degrees of freedom
## Multiple R-squared: 0.01212, Adjusted R-squared: -0.02316
## F-statistic: 0.3437 on 1 and 28 DF, p-value: 0.5624
# 近似直線も男女別に描きます
abline(fit.lm.2.m,col=4)
abline(fit.lm.2.f,col=2)
# 近似直線を見てみると、男女で傾きが違うように見えます。次の仮説を立てました。
# 「男では入学試験成績と卒業試験成績に相関がある。女では相関は無い。」言い換え
# ると、「入学試験成績と卒業試験成績の相関関係は性別によって修飾される（交互作
# 用がある）」という仮説になります。検証してみます。
fit.lm.3<- # 結果を格納するオブジェクト名は任意に決めます
lm( # 線形回帰分析ではlm関数を使います。
GradExam~EntranceExam*Sex, # 卒業試験成績を応答変数に、性別と入学試験成績を説明変
# 数にしますが、今回は*で結んでいます。*の意味は性別と
# 入学試験成績の交互作用項（EntranceExam:Sex）をモデル
# に含めなさいという意味です。
data=Dataset.1 # データセットの指定
)
summary(fit.lm.3)
##
## Call:
## lm(formula = GradExam ~ EntranceExam * Sex, data = Dataset.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.0698 -6.0168 0.4657 6.0183 22.4657
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63.8143 23.7604 2.686 0.0085 **
## EntranceExam 0.1269 0.2806 0.452 0.6523
## SexMale -24.7955 27.8710 -0.890 0.3758
## EntranceExam:SexMale 0.2374 0.3297 0.720 0.4731
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.151 on 98 degrees of freedom
## Multiple R-squared: 0.1117, Adjusted R-squared: 0.08446
## F-statistic: 4.106 on 3 and 98 DF, p-value: 0.008641
# 入学試験成績と性別の交互作用項もモデルに含めると、入学試験成績と卒業試験成績（P=0.652）、
# 男女差（P=0.376）、性別と入学試験成績の交互作用（P for interaction=0.473）のいずれも
# 有意差が無いことがわかります。有意差がないことは必ずしも「母集団に差がない」ことを示し
# ません。検出力が不足していても、モデルが間違っていても同じ結論を導きえます。
#
# モデルの違い、すなわち、卒業試験成績の予測に入学試験成績のみを用いたモデル（fit.lm.1）、
# 男女で卒業試験成績に差があると仮定したモデル（fit.lm.2）、入学試験成績と卒業試験成績の
# 関連は男女で異なると仮定したモデル（fit.lm.3）によって、仮定が異なれば当然ですが、結果
# が異なってくることに注意してください。
#
# 重回帰分析の理解は、全ての回帰分析と多変量解析の理解の基礎になります。この演習が理解の
# きっかけになることを期待します。