「なにかやりたい気持ち」改め、
「説明変数にノイズを含むモデルで遊んでみた」



この記事は Stan Advent Calendar 2016 の12月10日の記事です。
Stan Advent Calendar 2016 盛り上がってますね!!楽しいです!!

アヒル本 StanとRでベイズ統計モデリング の7.7「説明変数にノイズを含む」のモデルで遊んでみた、という内容です。モデルについての誤解がある気がしています。お気づきの点がございましたら、ご指摘いただければ幸いです。
 
 

説明変数にノイズを含むモデル


このモデルは説明変数をX、目的変数をYとしたときに、以下のようになります。
\begin{eqnarray}
X[n] & \sim & Normal(X_{true}[n], 1.0)\\
Y[n] & \sim & Normal(\alpha+\beta X_{true} [n], \sigma)
\end{eqnarray}


はじめの式は「説明変数が平均を真の値(Xtrue) 、標準偏差を1.0とした正規分布から発生している」、つまり「データとして得られた説明変数には標準偏差1.0のノイズが含まれている」ということを表現しています。ここではノイズについて1.0という固定値を与えていますが、もちろん変更できます。
二つ目の式では、説明変数の真の値で目的変数を回帰しています。アルファが切片、ベータが回帰係数、シグマが誤差標準偏差を示します。


構成概念をあつかう心理学では、説明変数に測定誤差などのノイズが含まれることが多いです。そこで複数項目で測定を行い、構造方程式モデリングで誤差と分離して目的変数への影響を検討しているわけなのですが、既存尺度を利用してデータを得た際に因子分析とかを各研究がそれぞれやり直しているような状況です。先行研究の因子分析の結果からは「どの項目にどの程度因子が寄与していて、ノイズがどのくらい含まれるのか」ということがわかるわけですから、そういった情報を利用して「少ない項目数、少ないサンプルサイズでも妥当な推定をする」ということができるといいなぁと思っていました。因子分析は潜在変数が説明変数、観測変数が目的変数となっている回帰分析と捉えることができます。はじめの式は「研究者が見たい真の値に誤差がのっかった状態で説明変数が生成されている」ということを意味していますから、アヒル本でこのモデルをみたときには興奮して、鼻血がでるんじゃないかと思いました。「そんなうまい話、、」というあれもあったので、鼻血はでなかったわけなのですが。


 

仮想データを発生させて遊んでみる


ということでデータの発生から。先に記述したモデルどおりです。

これをlmで分析すると下の結果。真の値 (t) を使って回帰分析をすると、(当たり前ですが) 切片も回帰係数も設定した値 (それぞれ3と5) に近い値になります。説明変数を使って分析をすると、切片は3.04、回帰係数は2.33となりました。

 
 
はたして、説明変数にノイズが含まれていることを仮定したモデルだと、真の値で回帰したときの値に近い推定ができるのでしょうか?発生させた仮想データに対して、説明変数にノイズが含まれていることを表現した以下のモデルで推定します。モデルのコードはそのまんまですね。ここでは説明変数のノイズの大きさに、データ発生の際に設定した値である1を固定値として与えていますが、実際に使うときには、先行研究の結果に基づいて分布を与えることになるのではないかと思います。

 

推定結果はこちら。見事にノイズを含む説明変数をつかってlmで分析したときと同じ値が得られており、回帰係数については設定した値である5とはほど遠い値が得られました。やはりそんなうまい話はないのか。。


 
 
ちょっと普通にX1で回帰分析するモデルで推定してみます。回帰分析における誤差の部分 (e) が、説明変数にノイズがあることを仮定するモデルよりも大きくなっています (ただし、lmで分析したときの残差の標準偏差は3.61なので、不当に大きいとか、そういうことではないです)。「説明変数にノイズが含まれている」という仮定は回帰係数の推定にではなく、誤差の推定に効いてくるようです。

 

よろしい、ならば3項目だ


変数が1つだけというのはさすがに無茶だったのかなと思い、変数を増やしてみました。データの発生は変数の数が増えただけで、先ほどと同様です。また、真の値 (t) を使って回帰分析をすると、やはり切片も回帰係数も設定した値 (それぞれ3と5) に近い値になります。3つの変数を平均した値を説明変数としてもちいると回帰係数は3.81となりました。


さて、このデータを以下のモデルで推定します。3項目の平均値をつかうモデルでも推定し、比較を行います。


結果はこちら。やはり、両モデルの違いは、切片や回帰係数の値ではなく、誤差の部分にあらわれています。


なお、このデータを普通にSEMで分析してみると、設定した値5に近い値が得られました。

 
 

まとめ


心理学のデータでは説明変数にノイズが含まれていることが多く、「ノイズを除去してより正確な変数間の推定ができれば最高だ!」と興奮したわけですが、回帰係数の推定に効くという感じではなさそうです (仮想データの発生の仕方が悪いのかなぁ)。このあたり、自分はまだまだモデルについての感覚が甘く、遊んでみないとわからないという部分でした。先行研究の因子分析の結果をつかって、「少ない項目数、少ないサンプルサイズでも妥当な推定をする」ということができるといいなぁと思ってます。

 
 

【追記】事前分布の設定をするとうまくいく!?


2016.12.12追記
関西学院大学の清水裕士先生に事前分布の設定をするとうまくいくというご助言をいただきました。
モデルにはパラメタの事前分布についての設定を追加しています。特に、パラメタtについては「平均mu、標準偏差sigmaに従う正規分布から発生している」という仮定をいれており、muとsigmaを新たにパラメタとして推定しています。

 
 

推定結果はこちら。データ発生の際に設定した値はa=3, b=5, e=0.5, mu=0, sigma=1なのですが、(特にtについての) 事前分布を設定しない場合に比べ、今回のモデルではそれぞれ近い値が得られています。ただし、今回は反復回数を10万、間引き数を5に設定していて、4万のmcmcサンプルがあるのですが、有効サンプル数 (n_eff) がとても少ないです。もともと説明変数にノイズが含まれていることを表現したモデルだと反復回数や間引き数を多くしないと収束していなかったのですが、事前分布の設定を追加すると、さらに収束のために反復回数や間引き数が必要な感じです。

One thought on “「なにかやりたい気持ち」改め、
「説明変数にノイズを含むモデルで遊んでみた」

  1. モデル式7-6のような課題を解くモデルを、変数内誤差モデルといいます。
    予測変数を3つにしていますが、予測変数が1つでもSEMのソフトで計算できます。
    具体例に関しましては、拙著「共分散構造分析、応用編(2000)、8章3節」
    ご参照ください。
    心理学の応用では、年齢の平均的な評価誤差が既知になるようなことは珍しく
    xの信頼性係数が別途推定されている場合がほとんどです。
    その場合の変数内誤差の扱いは、拙著「原因を探る統計学(1992)p.196-200」を
    ご参照ください。
    共分散構造をstanで解くことも可能ですので、それを以下に示します。
    要するに清水さんが示したコードのtを積分消去して周辺ベイズ推定しています。
    要するに共分散構造分析とは、
    そもそも潜在変数を積分消去した周辺ベイズ推定なのです。

    library(rstan)
    N<-300 #データ数
    a<-3 #切片
    b<-5 #傾き
    et<-10  #真の年齢の散らばり
    mu<-30 #写真の真の平均年齢
    ey<-1.5 #誤差SD
    ex<-2.5 #年齢の平均的な評価誤差(既知)
    set.seed(888)
    t<-rnorm(N,mu,et)
    x<-t + rnorm(N,0,ex)
    y<-a + b*t +rnorm(N,0,ey)
    xy<-cbind(x,y)

    #真の共分散
    print(et*et + ex*ex)
    print(b * et*et)
    print((b*b * et*et) + ey*ey)
    #データの共分散
    print(cov(xy))

    #真の平均
    print(mu)
    print(a+b*mu)
    #データの平均
    colMeans(xy)

    data<-list(N=N,xy=xy)
    fit<-stan(file="model7_6new.stan",data=data,
    pars = c("a","b","et","mu","ey","mea","cov"))
    print(fit)
    traceplot(fit, include =T)

    ーーーーーーーーーーーーーーーーーーーーーーーー"model7_6new.stan"
    data{
    int N;
    vector[2] xy[N];
    }

    transformed data{
    real ex;
    ex = 2.5; #年齢の平均的な評価誤差(既知)
    }

    parameters{
    real a;
    real b;
    real et;
    real mu;
    real ey;
    }

    transformed parameters{
    vector[2] mea;
    cov_matrix[2] cov;
    #共分散構造
    cov[1,1] =et*et + ex*ex;
    cov[1,2] =b * et*et;
    cov[2,1] =b * et*et;
    cov[2,2] =(b*b * et*et) + ey*ey;
    #平均構造
    mea[1] =mu;
    mea[2] =a+b*mu;
    }
    model{
    xy ~ multi_normal(mea,cov);
    }

    結果は以下で、aがズレテいますが、
    発生したデータの標本誤差の影響と思われます。

    post-warmup draws per chain=1000, total post-warmup draws=4000.

    mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
    a 4.57 0.06 2.40 -0.26 2.94 4.60 6.19 9.10 1705 1.00
    b 4.97 0.00 0.08 4.83 4.92 4.97 5.02 5.13 1727 1.00
    et 9.90 0.01 0.42 9.11 9.62 9.89 10.18 10.79 2675 1.00
    mu 30.34 0.01 0.56 29.24 29.96 30.35 30.72 31.43 2924 1.00
    ey 1.47 0.02 1.03 0.08 0.65 1.31 2.10 3.90 2624 1.00
    mea[1] 30.34 0.01 0.56 29.24 29.96 30.35 30.72 31.43 2924 1.00
    mea[2] 155.51 0.05 2.72 150.25 153.61 155.52 157.38 160.85 3017 1.00
    cov[1,1] 104.53 0.16 8.46 89.27 98.82 103.99 109.79 122.72 2674 1.00
    cov[1,2] 488.69 0.74 40.30 416.63 461.04 486.51 513.97 575.76 2973 1.00
    cov[2,1] 488.69 0.74 40.30 416.63 461.04 486.51 513.97 575.76 2973 1.00
    cov[2,2] 2433.83 3.47 198.00 2079.21 2295.87 2426.23 2558.07 2848.67 3247 1.00
    lp__ -1723.38 0.04 1.58 -1727.11 -1724.25 -1723.09 -1722.19 -1721.25 1378 1.01

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です