RStanで重回帰を試してみる - Gakushukun1’s diary

目的

e-statの統計データから, rstanによる重回帰を用いて年代ごとの趣味の関係を調べる.

今回はRStanを使って重回帰分析を試してみる. 対象としたのは, e-statの中の,
社会生活基本調査平成28年社会生活基本調査調査票Ａに基づく結果生活行動に関する結果生活行動編（全国）趣味・娯楽である. このデータでは, おおまかな世代ごとの趣味の傾向について調査の結果が掲載されている. 例として, 10〜14歳の世代のデータでは, 調査を行ったサンプル9351人のうち, 趣味にスポーツ観覧が含まれると答えた人が1804人, 美術鑑賞が含まれていると答えた人が854人ということが分かる(なお, これらは重複解答を含んでいる).

本エントリでは, 複数ある趣味のカテゴリの項目のうち, スポーツ観覧, 園芸, とテレビ・パソコンゲームの関係を調べてみることとする. 使用するモデルは次の通り.

$Z = a + bX + cY + \epsilon \\$

このモデルを用いて, $X$ をスポーツ観覧が趣味の人の割合, $Y$ を園芸が趣味の人の割合, $Z$ をテレビ・パソコンゲームが趣味の人の割合として回帰を行う. ここでいう割合とは, 調査を行った各世代ごとのサンプルサイズに対して, それぞれの趣味を持っていると回答した人の割合である. なお, 回帰に用いるデータは10〜14歳の世代, 15～19歳の世代, …, 85歳以上の世代の16個である.

重回帰分析に用いたコードを次に示す.

library(rstan)
library(maptools)

d <- read.csv("FEH_00200533_190529214744.csv", fileEncoding="shift-jis", stringsAsFactors=FALSE)
age_str <- strsplit(unlist(d[2:17,6]), "_")
age <- c()
for(i in 1:length(age_str)) {
    age <- append(age, age_str[[i]][2])
}
sample_size <- as.numeric(unlist(d[2:17,8]))
sport <- as.numeric(unlist(d[2:17,14])) / sample_size
garden <- as.numeric(unlist(d[2:17,54])) / sample_size
game <- as.numeric(unlist(d[2:17,74])) / sample_size
df <- data.frame(age, sport, garden, game)

data <- list(N=nrow(df), X=df$sport, Y=df$garden, Z=df$game)
fit <- stan(file="jukaiki.stan", data=data, seed=123)

ms <- extract(fit)
x <- seq(0.0, 0.3, length=50)
y <- seq(0.0, 0.3, length=50)
func <- function(x, y) {mean(ms$a) + mean(ms$b) * x + mean(ms$c) * y}
z <- outer(x, y, func)
png("hobby.png", width=2000, height=2000, res=300 )
par(family = "HiraKakuProN-W3")
contour(x, y, z)
par(new=T)
plot(df$sport, df$garden, xlab="watching_sports", ylab="gardening", xlim=c(0.0,0.3), ylim=c(0.0,0.3))
pointLabel(x=df$sport, y=df$garden, labels=df$age)
dev.off()

Stan "jukaiki.stan"

data {
    int N;
    real X[N];
    real Y[N];
    real Z[N];
}

parameters {
    real a;
    real b;
    real c;
    real<lower=0> sigma;
}

model {
    for (n in 1:N) {
        Z[n] ~ normal(a + b * X[n] + c * Y[n], sigma);
    }
}

generated quantities {
    vector[N] log_lik;

    for (n in 1:N) {
        log_lik[n] = normal_lpdf(Z[n] | a + b * X[n] + c * Y[n], sigma);
    }
}

各パラメータのサンプリング結果の平均は次のようになった.

a	b	c	sigma
0.07	2.08	-0.78	0.04

これらの値を元のモデルに当てはめ, テレビ・パソコンゲームを等高線でプロットした図は次のようになる. なお, 図中の点は分析に用いた元のデータをプロットしたものである. サンプルサイズが小さいものの, 良好に推測できていると考えられる.
f:id:gakushukun1:20190531010834p:plain