『Stataによる計量政治学』（浅野正彦, 矢内勇生)：Rによる分析
Yuki Yanai
October 10, 2014 (rev. June 6, 2015)
第13章 回帰分析の応用
13.1 ダミー変数の利用
13.1.1 ダミー変数
準備として、作業ディレクトリを指定する。 作業ディレクトリの指定にはsetwd()関数を利用する。 ディレクトリへのパス（引用符の中身）は、各自で設定する。
setwd("~/projects/Stata-R/")この章で使うパッケージをあらかじめ読み込んでおく。
library("foreign")
library("dplyr")##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary("ggplot2")
## Windows ユーザは次の行をコメントアウト
theme_set(theme_gray(base_size=12, base_family="HiraKakuProN-W3"))再び、2009年の衆院選における自民党候補の得票率を選挙費用と当選回数で説明する回帰モデルを例として使う。
まず、分析用のデータ (hr96-09.dta)を読み込む （第5章のRによる分析で作った hr96-09labeled.dat を使ってもよい）。
HR <- read.dta("data/hr96-09.dta")
## hr96-09labeled.dat を使う場合は、
# load("data/hr96-09labeled.dat")
## データが適切に読み込まれたどうか確認したいなら、
# head(HR)では、ダミー変数を作ろう。 ここでは、データセットに含まれる status というカテゴリー変数からダミー変数を作ってみよう。 まず、status 変数がどんな変数か確認する。
with(HR, table(status))## status
## challenger incumbent moto
## 3411 1908 295ここからわかる通り、status という変数は、「新人 (challenger)」「現職 (incumbent)」「元職 (moto)」を区別するカテゴリー変数である。
このカテゴリー変数から、新人候補を表すダミー変数 new と、元職をダミー示すダミー変数 old を作ろう。 ダミー変数は、as.numeric() 関数で作ることができる。 関数の引数には、ダミー変数の値が1になる条件を与える。
HR <- mutate(HR,
new = as.numeric(status == "challenger"),
old = as.numeric(status == "moto"))ここで、ダミー変数とその元となったカテゴリー変数を並べてみよう。
with(HR, cbind(status, new, old))上のコマンドを実行すると、3つの変数が並べて表示される（長くなるのでここでは省略する）。 3つの変数を比べると、newが1でoldがOのときはstatusが1、newがOでoldが1のときはstatusが3、newもold も0のときはstatusが2になることがわかる。 このように、カテゴリ一変数はその変数がもつカテゴリーの数より1つ少ない数のダミー変数で表現し直すことができる。
カテゴリー数が多いときは、for ループを利用してダミー変数を作る。 たとえば、party 変数の各カテゴリに対するダミー変数は次のようにして作ることができる。
lvs <- levels(HR$party) ## party変数のカテゴリ
n.party <- length(lvs) ## party変数のカテゴリ数
k <- dim(HR)[2] ## HRデータに含まれる変数の数
for (i in 1:n.party) { ## partyのカテゴリ数の分だけループする
HR[, k+i] <- as.numeric(HR$party == lvs[i]) ## i番目のカテゴリのダミー
names(HR)[k+i] <- paste("pty",lvs[i], sep=".") ## 変数名をつける
}これで、35個のダミー変数がHRデータに追加された。 次のコマンドで新しい変数を確認しておこう。
names(HR)## [1] "year" "ku" "kun" "party" "name"
## [6] "age" "status" "nocand" "wl" "rank"
## [11] "previous" "vote" "voteshare" "eligible" "turnout"
## [16] "exp" "new" "old" "pty.msz" "pty.JCP"
## [21] "pty.LDP" "pty.CGP" "pty.oki" "pty.tai" "pty.saki"
## [26] "pty.NFP" "pty.DPJ" "pty.SDP" "pty.LF" "pty.NJSP"
## [31] "pty.DRF" "pty.kobe" "pty.nii" "pty.sei" "pty.JNFP"
## [36] "pty.bunka" "pty.green" "pty.LP" "pty.RC" "pty.muk"
## [41] "pty.CP" "pty.NCP" "pty.ND" "pty.son" "pty.sek"
## [46] "pty.NP" "pty.NNP" "pty.NPJ" "pty.NPD" "pty.minna"
## [51] "pty.R" "pty.H" "pty.sho"## 新しいダミー変数の最初の6つの値
select(HR, pty.CGP:pty.sho) %>% head()## pty.CGP pty.oki pty.tai pty.saki pty.NFP pty.DPJ pty.SDP pty.LF pty.NJSP
## 1 0 0 0 0 1 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 1 0 0 0
## 4 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0
## pty.DRF pty.kobe pty.nii pty.sei pty.JNFP pty.bunka pty.green pty.LP
## 1 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 1 0 0
## 6 0 0 0 0 0 0 0 0
## pty.RC pty.muk pty.CP pty.NCP pty.ND pty.son pty.sek pty.NP pty.NNP
## 1 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 1 0
## pty.NPJ pty.NPD pty.minna pty.R pty.H pty.sho
## 1 0 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 0
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 0 0 0 0## 新しいダミー変数の最後の6つの値
select(HR, pty.CGP:pty.sho) %>% tail()## pty.CGP pty.oki pty.tai pty.saki pty.NFP pty.DPJ pty.SDP pty.LF
## 5609 0 0 0 0 0 0 0 0
## 5610 0 0 0 0 0 0 0 0
## 5611 0 0 0 0 0 0 0 0
## 5612 0 0 0 0 0 1 0 0
## 5613 0 0 0 0 0 0 0 0
## 5614 0 0 0 0 0 0 0 0
## pty.NJSP pty.DRF pty.kobe pty.nii pty.sei pty.JNFP pty.bunka
## 5609 0 0 0 0 0 0 0
## 5610 0 0 0 0 0 0 0
## 5611 0 0 0 0 0 0 0
## 5612 0 0 0 0 0 0 0
## 5613 0 0 0 0 0 0 0
## 5614 0 0 0 0 0 0 0
## pty.green pty.LP pty.RC pty.muk pty.CP pty.NCP pty.ND pty.son pty.sek
## 5609 0 0 0 0 0 0 0 0 0
## 5610 0 0 0 0 0 0 0 0 0
## 5611 0 0 0 0 0 0 0 0 0
## 5612 0 0 0 0 0 0 0 0 0
## 5613 0 0 0 0 0 0 0 0 0
## 5614 0 0 0 0 0 0 0 0 0
## pty.NP pty.NNP pty.NPJ pty.NPD pty.minna pty.R pty.H pty.sho
## 5609 0 0 0 0 0 0 0 0
## 5610 0 0 0 0 0 0 0 0
## 5611 0 0 0 0 0 0 1 0
## 5612 0 0 0 0 0 0 0 0
## 5613 0 0 0 0 0 0 0 0
## 5614 0 0 0 0 0 0 1 013.1.2 ダミー変数を使った回帰分析
ダミー変数を使って切片を変える
教科書 p.227 の例13-1 を使って説明する。
まず、衆院選データセットHRから、2009年のデータだけ取り出そう。
HR2009 <- filter(HR, year == 2009)民主党候補に興味があるので、民主党候補を表すダミー変数dpjをつくる。
HR2009 <- mutate(HR2009,
dpj = as.numeric(party == "DPJ"),
dpj = factor(dpj, labels = c("others", "DPJ")))このダミー変数を回帰分析で利用する。 教科書と同じように、得票率 (voteshare) を応答変数、選挙費用 (exp、ただし、100万で割って使う) と民主党ダミー (dpj) を説明変数とする回帰分析を行う。
fit.dpj <- lm(voteshare ~ I(exp / (10^6)) + dpj, data = HR2009)
summary(fit.dpj)##
## Call:
## lm(formula = voteshare ~ I(exp/(10^6)) + dpj, data = HR2009)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39.671 -6.712 -3.042 5.903 63.473
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.89321 0.53949 7.217 9.81e-13
## I(exp/(10^6)) 2.62981 0.06805 38.646 < 2e-16
## dpjDPJ 27.45274 0.79766 34.417 < 2e-16
##
## Residual standard error: 11.19 on 1121 degrees of freedom
## (15 observations deleted due to missingness)
## Multiple R-squared: 0.7463, Adjusted R-squared: 0.7458
## F-statistic: 1649 on 2 and 1121 DF, p-value: < 2.2e-16このように、表13.2 (p.228) と同じ結果が得られる（結果の解釈については教科書を参照）。
この結果を図示してみよう。
## 民主党候補を赤色、それ以外の候補を青色でプロットする
plot(voteshare ~ I(exp / (10^6)), data = HR2009,
main = "2009年衆院選の得票率：選挙費用と民主党ダミーで説明する",
ylab = "得票率 (%)", xlab = "選挙費用（100万円）",
col = ifelse(HR2009$dpj == "DPJ", "red", "blue"), # col で色を指定
pch = ifelse(HR2009$dpj == "DPJ", 5, 16)) # pch で点の形を指定
## 回帰直線を上書きする
## 直線を引く範囲を求める（データの範囲外に直線を引かない！）
xrng.dpj <- HR2009 %>%
filter(dpj == "DPJ", !is.na(exp)) %>%
with(c(min(exp), max(exp)) / 10^6)
xrng.oth <- HR2009 %>%
filter(dpj == "others", !is.na(exp)) %>%
with(c(min(exp), max(exp)) / 10^6)
## 民主党候補の回帰直線を赤色で描く
lines(xrng.dpj, coef(fit.dpj)[1] + coef(fit.dpj)[3] + coef(fit.dpj)[2] * xrng.dpj,
col = "red", lwd = 2)
## その他の候補の回帰直線を青色で描く
lines(xrng.oth, coef(fit.dpj)[1] + coef(fit.dpj)[2] * xrng.oth,
col = "blue", lwd = 2)
## 図の左上のスペースに凡例 (legend) をつける
legend("topleft", title="所属政党", legend = c("民主党", "その他"),
col = c("red", "blue"), pch = c(5, 16))
これで、図13.1 (p.229) と同じような図ができた（白黒で印刷する場合に備えて点の形 [pch] も変えている）。
次に、民主党に所属していたどうかだけでなく、所属政党ごとに投票率が異なるかどうか検討しよう（p.230, 例13-2）。 Rでは、
fit.party <- lm(voteshare ~ I(exp / (10^6)) + party, data = HR2009)
summary(fit.party)##
## Call:
## lm(formula = voteshare ~ I(exp/(10^6)) + party, data = HR2009)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.686 -3.052 -0.177 2.209 70.835
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.71789 1.08838 9.848 < 2e-16
## I(exp/(10^6)) 0.72305 0.07736 9.346 < 2e-16
## partyJCP -4.28399 1.23136 -3.479 0.000523
## partyLDP 21.65672 1.24865 17.344 < 2e-16
## partyCGP 20.62830 3.47081 5.943 3.73e-09
## partyDPJ 35.50573 1.16688 30.428 < 2e-16
## partySDP 7.25437 1.76503 4.110 4.25e-05
## partyNNP 21.55124 2.86390 7.525 1.08e-13
## partyNPJ 26.82623 5.73505 4.678 3.26e-06
## partyminna 6.06706 2.36967 2.560 0.010589
## partyR -4.37604 8.05205 -0.543 0.586916
## partyH -10.38753 1.15236 -9.014 < 2e-16
##
## Residual standard error: 7.978 on 1112 degrees of freedom
## (15 observations deleted due to missingness)
## Multiple R-squared: 0.8721, Adjusted R-squared: 0.8708
## F-statistic: 689.3 on 11 and 1112 DF, p-value: < 2.2e-16これで、表13.3 (p.232) と同じ結果が得られた（表の読み方については教科書を参照）。
（注：ここで利用しているカテゴリ変数 party のクラスはfactor である。あらかじめfactor として用意されていないカテゴリ変数を利用するときは、as.factor() 関数を用い、 lm(votershare ~ I(exp / (10^6) + as.factor(party), data = HR2009) のようにする。変数のクラスは、class(変数名) で確認できる。）
この分析では、基準カテゴリ（参照カテゴリ）は無所属の議員である。 基準となるカテゴリは、利用するカテゴリ変数によって決まる。 分析で利用したカテゴリ変数の中身を確かめてみると、
levels(HR2009$party)## [1] "msz" "JCP" "LDP" "CGP" "oki" "tai" "saki" "NFP"
## [9] "DPJ" "SDP" "LF" "NJSP" "DRF" "kobe" "nii" "sei"
## [17] "JNFP" "bunka" "green" "LP" "RC" "muk" "CP" "NCP"
## [25] "ND" "son" "sek" "NP" "NNP" "NPJ" "NPD" "minna"
## [33] "R" "H" "sho"となっている。 この変数は、msz（無所属）というカテゴリが第1のカテゴリになるように作られている。 そして、この第1のカテゴリが回帰分析の基準カテゴリになっている。
基準となるカテゴリを変えたいときは、新たなカテゴリ変数を作ればよい。 自民党 (LDP) を基準カテゴリにした分析を行ってみよう。
まず、自民党カテゴリが第1カテゴリになるようなカテゴリ変数（factor クラスの変数）を作る。 既に存在するparty 変数に対してfactor() 関数を適用し、カテゴリの順番を指定してやればよい。 ここでは、初めに自民党 (LDP) を指定し、その後の順番はparty と同じ（ただし、元々3番目にあったLDPを飛ばすため、“-3” とする）にする。
HR2009 <- mutate(HR2009,
party.LDPbase = factor(party, c("LDP", levels(party)[-3])))新しい変数のカテゴリを確認してみると、
levels(HR2009$party.LDPbase)## [1] "LDP" "msz" "JCP" "CGP" "oki" "tai" "saki" "NFP"
## [9] "DPJ" "SDP" "LF" "NJSP" "DRF" "kobe" "nii" "sei"
## [17] "JNFP" "bunka" "green" "LP" "RC" "muk" "CP" "NCP"
## [25] "ND" "son" "sek" "NP" "NNP" "NPJ" "NPD" "minna"
## [33] "R" "H" "sho"というように、希望通りの順番の変数ができている。
この変数を使って回帰分析を行ってみよう。
fit.party2 <- lm(voteshare ~ I(exp / (10^6)) + party.LDPbase, data = HR2009)
summary(fit.party2)##
## Call:
## lm(formula = voteshare ~ I(exp/(10^6)) + party.LDPbase, data = HR2009)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.686 -3.052 -0.177 2.209 70.835
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.37461 0.99694 32.474 < 2e-16
## I(exp/(10^6)) 0.72305 0.07736 9.346 < 2e-16
## party.LDPbasemsz -21.65672 1.24865 -17.344 < 2e-16
## party.LDPbaseJCP -25.94071 1.07538 -24.122 < 2e-16
## party.LDPbaseCGP -1.02842 3.29158 -0.312 0.75476
## party.LDPbaseDPJ 13.84901 0.73029 18.964 < 2e-16
## party.LDPbaseSDP -14.40235 1.58485 -9.088 < 2e-16
## party.LDPbaseNNP -0.10548 2.71313 -0.039 0.96900
## party.LDPbaseNPJ 5.16951 5.67768 0.910 0.36276
## party.LDPbaseminna -15.58966 2.22610 -7.003 4.32e-12
## party.LDPbaseR -26.03276 7.99327 -3.257 0.00116
## party.LDPbaseH -32.04425 0.99908 -32.074 < 2e-16
##
## Residual standard error: 7.978 on 1112 degrees of freedom
## (15 observations deleted due to missingness)
## Multiple R-squared: 0.8721, Adjusted R-squared: 0.8708
## F-statistic: 689.3 on 11 and 1112 DF, p-value: < 2.2e-16これで、表13.4 (p.233) と同じ結果が得られた。
ダミー変数を使って切片と傾きを変える
続いて、選挙費用が得票率に与える影響が、民主党候補とその他の候補で異なるかどうか検証する（p.234, 例13-3）。 つまり、選挙費用と民主党ダミーの交差項を回帰モデルに導入する。
fit.dpj.int <- lm(voteshare ~ I(exp / (10^6)) * dpj, data = HR2009)
summary(fit.dpj.int)##
## Call:
## lm(formula = voteshare ~ I(exp/(10^6)) * dpj, data = HR2009)
##
## Residuals:
## Min 1Q Median 3Q Max
## -41.534 -5.755 -2.464 4.383 62.031
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.29096 0.52339 4.377 1.31e-05
## I(exp/(10^6)) 2.91640 0.06807 42.845 < 2e-16
## dpjDPJ 45.77138 1.67183 27.378 < 2e-16
## I(exp/(10^6)):dpjDPJ -2.42900 0.19817 -12.257 < 2e-16
##
## Residual standard error: 10.51 on 1120 degrees of freedom
## (15 observations deleted due to missingness)
## Multiple R-squared: 0.7763, Adjusted R-squared: 0.7757
## F-statistic: 1296 on 3 and 1120 DF, p-value: < 2.2e-16これで、表13.5 (p.236) と同じ結果が得られた。 このように、交差させる（相互作用を考える）変数同士を+ の代わりに* で繋ぐと、切片と傾きの両者をダミー変数の値によって変えるようなモデルを使った分析ができる。
切片と傾きの両者をダミー変数を使って変えるようなモデルを図示するときは、ggplot2 を使うのが便利である。
fig13.2 <- ggplot(HR2009, aes(I(exp /( 10^6)), voteshare, color = dpj, shape = dpj))
fig13.2 + geom_point() + geom_smooth(method = "lm", se = FALSE) +
xlab("選挙費用（100万円）") + ylab("得票率 (%)") +
scale_color_discrete(name = "政党", labels = c("その他", "民主党")) +
scale_shape_discrete(name = "政党", labels = c("その他", "民主党")) +
guides(color = guide_legend(reverse = TRUE), shape = guide_legend(reverse = TRUE))
これで、図13.2 (p.237) と同種の図ができた（白黒で印刷する場合に備えて点の形 [shape] も変えている）。
13.2 変数変換
例13-3「父親の身長が高いほど、その子どもの身長は高くなるか」(p.240) について、教科書と同じ架空のデータ（height.dta）を使って考える。
まず、データを読み込んで、データセットの概要を把握しよう。
HT <- read.dta("data/height.dta")
head(HT)## ht female father metfather infather
## 1 145.4 1 173.3 1.733 68.22835
## 2 158.8 1 165.6 1.656 65.19685
## 3 171.6 0 184.0 1.840 72.44095
## 4 186.0 0 175.1 1.751 68.93701
## 5 188.6 0 165.5 1.655 65.15748
## 6 157.1 1 170.6 1.706 67.16536str(HT)## 'data.frame': 100 obs. of 5 variables:
## $ ht : num 145 159 172 186 189 ...
## $ female : int 1 1 0 0 0 1 1 0 0 0 ...
## $ father : num 173 166 184 175 166 ...
## $ metfather: num 1.73 1.66 1.84 1.75 1.65 ...
## $ infather : num 68.2 65.2 72.4 68.9 65.2 ...
## - attr(*, "datalabel")= chr ""
## - attr(*, "time.stamp")= chr "24 Dec 2012 14:53"
## - attr(*, "formats")= chr "%9.0g" "%8.0g" "%9.0g" "%9.0g" ...
## - attr(*, "types")= int 254 251 254 254 254
## - attr(*, "val.labels")= chr "" "" "" "" ...
## - attr(*, "var.labels")= chr "height in cm" "dummy indicating female" "father's height in cm" "father's height in meter" ...
## - attr(*, "version")= int 12summary(HT)## ht female father metfather
## Min. :142.2 Min. :0.00 Min. :153.6 Min. :1.536
## 1st Qu.:157.3 1st Qu.:0.00 1st Qu.:165.6 1st Qu.:1.656
## Median :164.1 Median :0.00 Median :168.8 Median :1.688
## Mean :164.6 Mean :0.48 Mean :169.1 Mean :1.691
## 3rd Qu.:171.4 3rd Qu.:1.00 3rd Qu.:172.3 3rd Qu.:1.723
## Max. :188.6 Max. :1.00 Max. :184.0 Max. :1.840
## infather
## Min. :60.47
## 1st Qu.:65.19
## Median :66.46
## Mean :66.59
## 3rd Qu.:67.82
## Max. :72.44データの特徴がつかめたところで、説明変数を父親の身長 (father) と女性ダミー (female)、応答変数を子の身長 (ht) とする回帰分析を行う。
fit.ht1 <- lm(ht ~ father + female, data = HT)
summary(fit.ht1)##
## Call:
## lm(formula = ht ~ father + female, data = HT)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.9745 -5.1455 -0.5746 5.3316 23.6153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.3505 27.5524 3.933 0.000158
## father 0.3631 0.1619 2.243 0.027187
## female -10.7016 1.9327 -5.537 2.63e-07
##
## Residual standard error: 9.552 on 97 degrees of freedom
## Multiple R-squared: 0.2929, Adjusted R-squared: 0.2783
## F-statistic: 20.09 on 2 and 97 DF, p-value: 5.007e-08これで、表13.6 (p.241) と同じ結果が得られた。
13.2.1 線形変換
教科書と同じように父親の身長を
- cmで測定した father
- meter で測定した metfather
- inch で測定した infather
とうい3つの説明変数を使って回帰分析をしてみよう（1は既に推定済みなので、上で求めたものを使う）。
fit.ht2 <- lm(ht ~ metfather + female, data = HT)
fit.ht3 <- lm(ht ~ infather + female, data = HT)
summary(fit.ht1) ## 父親の身長の測定単位 = cm##
## Call:
## lm(formula = ht ~ father + female, data = HT)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.9745 -5.1455 -0.5746 5.3316 23.6153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.3505 27.5524 3.933 0.000158
## father 0.3631 0.1619 2.243 0.027187
## female -10.7016 1.9327 -5.537 2.63e-07
##
## Residual standard error: 9.552 on 97 degrees of freedom
## Multiple R-squared: 0.2929, Adjusted R-squared: 0.2783
## F-statistic: 20.09 on 2 and 97 DF, p-value: 5.007e-08summary(fit.ht2) ## 父親の身長の測定単位 = m##
## Call:
## lm(formula = ht ~ metfather + female, data = HT)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.9745 -5.1455 -0.5746 5.3316 23.6153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.350 27.552 3.933 0.000158
## metfather 36.313 16.191 2.243 0.027187
## female -10.702 1.933 -5.537 2.63e-07
##
## Residual standard error: 9.552 on 97 degrees of freedom
## Multiple R-squared: 0.2929, Adjusted R-squared: 0.2783
## F-statistic: 20.09 on 2 and 97 DF, p-value: 5.007e-08summary(fit.ht3) ## 父親の身長の測定単位 = inch##
## Call:
## lm(formula = ht ~ infather + female, data = HT)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.9745 -5.1455 -0.5746 5.3316 23.6153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.3505 27.5524 3.933 0.000158
## infather 0.9224 0.4113 2.243 0.027187
## female -10.7016 1.9327 -5.537 2.63e-07
##
## Residual standard error: 9.552 on 97 degrees of freedom
## Multiple R-squared: 0.2929, Adjusted R-squared: 0.2783
## F-statistic: 20.09 on 2 and 97 DF, p-value: 5.007e-08教科書で説明されているとおり、どの変数を使っても実質的に同じ結果が得られるが、推定値そのものは異なるので解釈のし易さに違いが出る。
Stata で分析する場合は上の例のようにあらかじめ変換された変数を用意して分析する必要があるが、Rの場合には分析するときに線形変換すればよい。 仮に、メートル単位で測ったmetfather しかデータセットに含まれていないとき、センチメートル単に直して分析したいときは、I() を使って次のようなコードを書けばよい。
fit.ht4 <- lm(ht ~ I(100 * metfather) + female, data = HT)
summary(fit.ht4)##
## Call:
## lm(formula = ht ~ I(100 * metfather) + female, data = HT)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.9745 -5.1455 -0.5746 5.3316 23.6153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.3505 27.5524 3.933 0.000158
## I(100 * metfather) 0.3631 0.1619 2.243 0.027187
## female -10.7016 1.9327 -5.537 2.63e-07
##
## Residual standard error: 9.552 on 97 degrees of freedom
## Multiple R-squared: 0.2929, Adjusted R-squared: 0.2783
## F-statistic: 20.09 on 2 and 97 DF, p-value: 5.007e-08これは、father変数を使って分析したfit.ht1 とまったく同じである。
13.2.2 中心化
標本平均を使って変数の中心化 (centering) を行う。 father を中心化したcfather と、female を中心化したcfemale という変数をデータセットに加えてみよう。
HT <- mutate(HT,
cfather = father - mean(father),
cfemale = female - mean(female))
mean(HT$cfather)## [1] 9.095019e-15mean(HT$cfemale)## [1] 1.775056e-17中心化された変数の平均値は実質的には0になっている。
これらの中心化された変数を使って回帰分析を行ってみよう。
fit.ht.ctr <- lm(ht ~ cfather + cfemale, data = HT)
summary(fit.ht.ctr)##
## Call:
## lm(formula = ht ~ cfather + cfemale, data = HT)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.9745 -5.1455 -0.5746 5.3316 23.6153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 164.6320 0.9552 172.353 < 2e-16
## cfather 0.3631 0.1619 2.243 0.0272
## cfemale -10.7016 1.9327 -5.537 2.63e-07
##
## Residual standard error: 9.552 on 97 degrees of freedom
## Multiple R-squared: 0.2929, Adjusted R-squared: 0.2783
## F-statistic: 20.09 on 2 and 97 DF, p-value: 5.007e-08これで、表13.7 (p.246) と同じ結果が得られた。