# ←この記号の後は次の改行までコメントになります。
# 今回の目標は、重回帰分析の点推定値、95%信頼区間、P-valueをsummary objectから
# 取り出して作表することです。
# 医学部学生の性別（Sex, 二項の名義変数）、入学成績（EntranceExam, 連続変数）、卒業
# 試験成績（GradExam, 連続変数）、医師国家試験の結果（NatnExam, 二項の名義変数）を
# 含む架空データを用います。
# ここから35行目まではサンプルデータセットの作成です。
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")
)
# まずは前回と同じ回帰分析を行います。
# 入学試験成績と性別から卒業試験成績を回帰分析する重回帰分析を用います。
# 再び回帰分析をしてみます。今回は多変量の重回帰分析ですが、文法はほぼ一緒です。
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
# Estimateが点推定値
# Std. Errorが標準誤差
# Pr(>|t|)がP-valueです。
# まず作業をスムースにするために、summary objectを保存します。
fit.lm.2.summary<-summary(fit.lm.2)
# summary objectの中身はlistです。中身を見るにはstr関数を使います。
str(fit.lm.2.summary)
## List of 11
## $ call : language lm(formula = GradExam ~ EntranceExam + Sex, data = Dataset.1)
## $ terms :Classes 'terms', 'formula' language GradExam ~ EntranceExam + Sex
## .. ..- attr(*, "variables")= language list(GradExam, EntranceExam, Sex)
## .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:3] "GradExam" "EntranceExam" "Sex"
## .. .. .. ..$ : chr [1:2] "EntranceExam" "Sex"
## .. ..- attr(*, "term.labels")= chr [1:2] "EntranceExam" "Sex"
## .. ..- attr(*, "order")= int [1:2] 1 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(GradExam, EntranceExam, Sex)
## .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "numeric" "factor"
## .. .. ..- attr(*, "names")= chr [1:3] "GradExam" "EntranceExam" "Sex"
## $ residuals : 'AsIs' Named num [1:102] -13.8164.... 2.110612.... -8.72415.... 8.110612.... 6.770345.... ...
## ..- attr(*, "names")= chr [1:102] "1" "2" "3" "4" ...
## $ coefficients : num [1:3, 1:4] 49.276 0.299 -4.763 12.503 0.147 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "(Intercept)" "EntranceExam" "SexMale"
## .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
## $ aliased : Named logi [1:3] FALSE FALSE FALSE
## ..- attr(*, "names")= chr [1:3] "(Intercept)" "EntranceExam" "SexMale"
## $ sigma : num 8.13
## $ df : int [1:3] 3 99 3
## $ r.squared : num 0.107
## $ adj.r.squared: num 0.0889
## $ fstatistic : Named num [1:3] 5.93 2 99
## ..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf"
## $ cov.unscaled : num [1:3, 1:3] 2.364854 -0.027592 -0.046746 -0.027592 0.000327 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "(Intercept)" "EntranceExam" "SexMale"
## .. ..$ : chr [1:3] "(Intercept)" "EntranceExam" "SexMale"
## - attr(*, "class")= chr "summary.lm"
# たくさんの出力がどっと出てきました。たとえば$callには、
fit.lm.2.summary$call
## lm(formula = GradExam ~ EntranceExam + Sex, data = Dataset.1)
# lm関数への入力が記載されています。
# これを全部理解する必要はとりあえずありません。まず必要なのは$coefficients
# だけです。見てみましょう。
fit.lm.2.summary$coef #他のリストの名前と区別がつけば、名前は$coefと短縮可
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.2763412 12.5034733 3.941012 0.0001513333
## EntranceExam 0.2988993 0.1469236 2.034387 0.0445867933
## SexMale -4.7630351 1.7683012 -2.693566 0.0083040945
# Estimateはcoefの1列目、Std. Errorは2列目、Pr(>|t|)は4列目に格納されてます。
# このうち必要なのは、EntranceExam（2行目）とSexMale（3行目）の値です。
fit.lm.2.summary$coef[c(2,3),1]
## EntranceExam SexMale
## 0.2988993 -4.7630351
fit.lm.2.summary$coef[c(2,3),2]
## EntranceExam SexMale
## 0.1469236 1.7683012
fit.lm.2.summary$coef[c(2,3),4]
## EntranceExam SexMale
## 0.044586793 0.008304094
# まずは、入学試験成績と卒業試験成績の関連を示す点推定値と95%信頼区間を
# 求めてみます。
# 95%信頼区間は点推定値 +/- 標準誤差*1.96です。
# 厳密には、正規分布の確率分布の2.5%と97.5%の値ですので、
# この正規分布の確率分布を与えるqnorm関数を用いて、
row.1<-
fit.lm.2.summary$coef[2,1]*c(1,1,1)+
# 点推定値を３つ並べる
fit.lm.2.summary$coef[2,2]*c(0,qnorm(0.025),qnorm(0.975))
# 後ろの2つに-1.96と1.96を標準誤差に掛けた値を加える。
row.1
## [1] 0.29889931 0.01093444 0.58686419
# 出てきた３つの値の最初が点推定値、後ろの２つが95%信頼区間の下限と上限です。
# これでは見にくいので、以前やったようにround関数とformat関数とpaste関数を用いて整形します。
c(
format(round(row.1[1],2),nsmall=2),
paste(
"[",
format(round(row.1[2],2),nsmall=2),
", ",
format(round(row.1[3],2),nsmall=2),
"]",
sep=""))
## [1] "0.30" "[0.01, 0.59]"
# もう一度コメント無しで繰り返して、今度はP-valueもくっつけましょう。
# また入学試験成績と性別の回帰分析を一度にやります。
row.1<-
c(fit.lm.2.summary$coef[2,1]*c(1,1,1)+
fit.lm.2.summary$coef[2,2]*c(0,qnorm(0.025),qnorm(0.975)),
fit.lm.2.summary$coef[2,4])
c(
format(round(row.1[1],2),nsmall=2),
paste(
"[",
format(round(row.1[2],2),nsmall=2),
", ",
format(round(row.1[3],2),nsmall=2),
"]",
sep=""),
format(round(row.1[4],3),nsmall=3))
## [1] "0.30" "[0.01, 0.59]" "0.045"
# 同じ計算を性別についても繰り返して、row.2を求め、
row.2<-
c(fit.lm.2.summary$coef[3,1]*c(1,1,1)+
fit.lm.2.summary$coef[3,2]*c(0,qnorm(0.025),qnorm(0.975)),
fit.lm.2.summary$coef[3,4])
# Tableを作ります。
Table<-
rbind(
c(
format(round(row.1[1],2),nsmall=2),
paste(
"[",
format(round(row.1[2],2),nsmall=2),
", ",
format(round(row.1[3],2),nsmall=2),
"]",
sep=""),
format(round(row.1[4],3),nsmall=3)),
c(
format(round(row.2[1],2),nsmall=2),
paste(
"[",
format(round(row.2[2],2),nsmall=2),
", ",
format(round(row.2[3],2),nsmall=2),
"]",
sep=""),
format(round(row.2[4],3),nsmall=3))
)
rownames(Table)<-c("Entrance examination per point","Male sex")
colnames(Table)<-c("Estimate","95%CI","P")
Table
## Estimate 95%CI P
## Entrance examination per point "0.30" "[0.01, 0.59]" "0.045"
## Male sex "-4.76" "[-8.23, -1.30]" "0.008"
# でもこのやりかたは、理解しておくと自在にTableを作ることができますが、
# ものすごく面倒くさいです。そこで、吉田和樹先生作のtableone::ShowRegTable
# を使って楽しましょう。
library(tableone)
print(ShowRegTable(fit.lm.2, # 用いた回帰分析の結果を使用
exp=F, # 線形回帰分析ではFを指定
ciFun=confint.default))
## coef [confint] p
## (Intercept) 49.28 [24.77, 73.78] <0.001
## EntranceExam 0.30 [0.01, 0.59] 0.045
## SexMale -4.76 [-8.23, -1.30] 0.008
## coef [confint] p
## (Intercept) "49.28 [24.77, 73.78]" "<0.001"
## EntranceExam " 0.30 [0.01, 0.59]" " 0.045"
## SexMale "-4.76 [-8.23, -1.30]" " 0.008"
# 95%信頼区間の計算方法
# たったこれだけ！