Hatena::ブログ(Diary)

アイアナ:データ分析や人工知能(AI)などの技術雑記

2009 | 11 |
2010 | 02 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2011 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2012 | 01 | 03 | 04 | 05 | 10 | 12 |
2013 | 01 | 02 | 04 |
2014 | 03 | 12 |
2016 | 03 |
2017 | 10 |

2011-05-23

小技メモ(Rの結果をHTMLで出力する)

今回の記事はここを参考にさせて頂きました。

http://d.hatena.ne.jp/dichika/20110520


Rは結果をレポートにするのが苦手で、ワードとかメモ帳に張り付けたりして編集したりするので結構な手間がかかります。

その点SPSSSASはレポートをそのままPDFとかHTMLに保存できるので楽だったりします。


実務レベルで統計ソフトを選ぶ場合は、レポート作成にかかる時間が馬鹿になりません。

その時間を節約するために、SASSPSSを解析ソフトとして選択するということは十分にあり得ました。


しかし、hwriterというパッケージを使えばRにもレポート利用の光が見えます!

これを使えば、テキスト・表・プロットもまとめてHTMLファイルに保存することができます。

library(hwriter)

#---準備
jpeg("test.jpg")
  plot(1:10)
dev.off()

a <- hwriteImage("test.jpg")
b <- hwrite("表です")
c <- hwrite(matrix(2, 5, 10))


#---ここからHTMLへ出力
p <- openPage("report.html", charset="shift-jis")
hwrite(a, p)
hwrite(b, p)
hwrite(c, p)
closePage(p)

これでreport.htmlがワーキングディレクトリデフォルトではマイドキュメンド)に保存されます。

デフォルトでは配置や表がアレですが、基本的にはHTMLなのでいろいろいじれば綺麗になりそうです。


javascriptcssも対応しているという噂も。。。

これを利用すれば、Rで解析→結果をwebページに直接アップということも出来そうです。


また改良出来次第、記事に致します。


追記

cssはopenPageのlink.cssオプションで指定できるようです!(@dichikaさんより)

2011-05-20

操作変数(instrumental variable, IV)での交絡調整

傾向スコアでの調整は、未測定の交絡要因は調整できないという問題がありました。

そこで操作変数(IV)という方法が提案されており、この方法では未測定の交絡要因も調整できると言われています。

それによって、よりランダム化試験に近い結果が出ると期待されます。


操作変数の定義は次のようになります。

  • A variable that is related to treatment but neither directly nor indirectly related to outcome, except through the effect of the treatment itself
  • 治療には関連しているが、結果には治療を通してでしか直接的にも間接的にも関連していない変数
    • 治療群Xを予測できる変数Zは、Xを通してでしか結果Yに関連しない(未知の交絡要因Uも介さない)

Rで実行するには、パッケージSEMを使います。

IVは計量経済(econometrics)の分野でも使われており、構造方程式モデルを利用してパラメータ推定を行います。

具体的には次のような2段階モデルを組みます。

  • X = α0 + α1 * Z + α2 * C + ε1
  • Y = β0 + β1 * X + β2 * C + ε2
    • X:治療、Y:結果、C: (複数の)測定済み交絡要因、Z:操作変数
    • αi, βi:係数
    • ε1, ε2:誤差(2変量正規分布を仮定することが多い)

プログラムの実行例はこのようになります。

library(sem)

data(Klein)
Klein$P.lag <- c(NA, Klein$P[-22])
Klein$X.lag <- c(NA, Klein$X[-22])
summary(tsls(C ~ P + P.lag + I(Wp + Wg), instruments=~1 + G + T + Wg + I(Year - 1931) + K.lag + P.lag + X.lag, data=Klein))

> 2SLS Estimates
>
>Model Formula: C ~ P + P.lag + I(Wp + Wg)
>
>Instruments: ~1 + G + T + Wg + I(Year - 1931) + K.lag + P.lag + X.lag
>
>Residuals:
>   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
> -1.890  -0.616  -0.246   0.000   0.885   2.000 
>
>            Estimate Std. Error t value  Pr(>|t|)
>(Intercept)  16.5548    1.46798 11.2772 2.587e-09
>P             0.0173    0.13120  0.1319 8.966e-01
>P.lag         0.2162    0.11922  1.8137 8.741e-02
>I(Wp + Wg)    0.8102    0.04474 18.1107 1.505e-12
>
>Residual standard error: 1.1357 on 17 degrees of freedom

tsls関数の「instruments=」の部分で操作変数のモデルを指定します。


この関数は結果変数が連続変数になっています。

2値結果変数の場合は直接tsls関数を直接利用できませんが、論文ベースでこのような提案がされています。

http://aje.oxfordjournals.org/content/169/3/273.abstract


ちなみに、未調整の線形モデルではこのような結果となります。

lm(C ~ P + P.lag + I(Wp + Wg), data = Klein)

>Coefficients:
>(Intercept)            P        P.lag   I(Wp + Wg)  
>   16.23660      0.19293      0.08988      0.79622  

調整することでPの効果が小さくなり、P.lagの効果が大きくなっているようです。

2011-05-16

Mongo DBを触ってみる

少し前からnoSQLのmongo DBが人気ということで少し導入してみました。

環境は64bit windowsです。


まずはダウンロード

http://www.mongodb.org/downloads


中身を任意の場所に解凍(自分はc:\mongo)。

モンゴが使うDBフォルダの作成。

  1. c:\data
  2. c:\data\db

cmdを実行してc:\mongo\binに移動。


mongodでモンゴサーバーが立ち上がる。

そのままの状態でhttp://localhost:28017/ブラウザするとこのようにモンゴサーバーを確認できる。

f:id:isseing333:20110516204217p:image


そのままcmdをもう一つ実行し、c:\mongo\binに移動してmongoと打つとめでたくmongo DBが動きます。

あとはここのチュートリアル通りやればモンゴの雰囲気が掴めますね〜

http://www.mongodb.org/pages/viewpage.action?pageId=5079135


サーバーはctrl+Cで終了することができます。

このチュートリアルのように使い方を分かりやすく訳された方々は素晴らしいと思います!

2011-05-11

観察データでの効果推定(傾向スコア、IPW、DR)

まずは教科書を紹介します。

Observational Studies (Springer Series in Statistics)

Observational Studies (Springer Series in Statistics)

Rosenbaum先生は傾向スコア(propensity score)を提案した方です。

この教科書に書いていあるのは傾向スコアについてだけで、IPWやDRは書いてありません。


日本語はこちらの星野先生の本。


IPWやDRのことも書いてあり、とても充実した内容になっています(今回とても参考になりました)。



観察研究はランダム化研究と違い、治療(曝露)対象者を選ぶことができません。

よって様々な交絡要因が混入している可能性があり、場合によっては調整する必要があります。

しかし調整できるのは、データとして収集している変数のみで、未測定の交絡要因は調整できません(ランダム化試験は理論上未測定の交絡要因も調整されています)。


交絡要因の調整方法は、次のようなものがあります。

  1. 傾向スコア
    • restriction(制限、調整)
    • stratification(層別)
    • matching(マッチング)
  2. IPW(Inverse Probability Weighting)
  3. DR推定量(ダブリーロバスト推定量、doubly robust)

Rのパッケージはあまり見当たりませんが、意外に簡単にRで推定できます。

データはMatchingパッケージに入っているlalondeを使います。



治療を割り付けられる「傾向」を使って調整する方法です。

まずロジスティック回帰で、ある対象者が治療が割り付けられる確率を推定します。

library(Matching)
data(lalonde)

#------傾向スコアの計算
Logit <- glm(treat ~ age + educ + black + hisp + married + nodegr + re74 + re75 + u74 + u75, 
	family = binomial, data = lalonde)
PS    <- Logit$fitted

次のように傾向スコアを使った調整(restriction)を行います。

#---restriction(制限、調整)
Restriction <- lm(lalonde$re78 ~ lalonde$treat + PS)
summary(Restriction)

これだけですので、傾向スコアでの調整は難しくありません。

ただ、データによっては効果が逆転してしまったり、標準誤差が大きくなってしまったりします。

またPSでは群間差は推定できるのですが、各群の期待値は推定できません。



  • IPW

IPWは傾向スコアPS)を利用して、「その治療群に割り付けられる確率」の逆数で重み付け推定します。

#------IPW(「その治療群に割り付けられる確率」の逆数で重み付けする)
Weight   <- lalonde$treat/PS + (1 - lalonde$treat)/(1 - PS)
IPweight <- lm(lalonde$re78 ~ lalonde$treat, weights = Weight)
summary(IPweight)

Y    <- lalonde$re78
(IPW1 <- sum( (Y/PS)[lalonde$treat == 1])/sum(1/PS[lalonde$treat == 1]))
(IPW0 <- sum( (Y/(1-PS))[lalonde$treat == 0])/sum(1/(1-PS)[lalonde$treat == 0]))


  • DR

ダブリーロバスト推定量はIPWをさらに拡張したもので、「2種のモデルのどちらかが間違っていてもバイアスのない推定ができる」推定方法です。

次の2モデルです。

  1. 割り付けに関するロジスティックモデル
  2. 結果変数に関する線形モデル

#------DR(ダブリーロバスト推定量)
n1      <- nrow(lalonde[lalonde$treat == 1, ])
n0      <- nrow(lalonde[lalonde$treat == 0, ])
Linear1 <- lm(re78 ~ age + educ + black + hisp + married + nodegr + re74 + re75 + u74 + u75, data = lalonde[lalonde$treat == 1, ])
Linear0 <- lm(re78 ~ age + educ + black + hisp + married + nodegr + re74 + re75 + u74 + u75, data = lalonde[lalonde$treat == 0, ])
Y1      <- lalonde$re78[lalonde$treat == 1]
Y0      <- lalonde$re78[lalonde$treat == 0]
PS1     <- PS[lalonde$treat == 1]
PS0     <- PS[lalonde$treat == 0]

(DR1 <- 1/n1 * sum(Y1 + ( (1 - PS1)/PS1) * (Y1 - Linear1$fitted)))
(DR0 <- 1/n0 * sum(Y0/(1 - PS0) + (1 - 1/(1 - PS0))*Linear0$fitted))


これらの推定量にどのような違いが出ているかをプロットします。

#------未調整の各群平均値
(Ttest <- t.test(lalonde$re78[lalonde$treat == 1], lalonde$re78[lalonde$treat == 0]))

#------差にエラーバーを付ける
#---Ttest
Ttestmean <- Ttest$estimate[1] - Ttest$estimate[2]
TtestCI <- Ttest$conf.int[1:2]

#---PS restriction
PsEst    <- summary(Restriction)$coefficients
Psmean   <- PsEst[2, 1]
Psse     <- PsEst[2, 2]
PsRestCI <- c(Psmean - 1.96*Psse, Psmean + 1.96*Psse)

#---IPW
IpwEst   <- summary(IPweight)$coefficients
Ipwmean  <- IpwEst[2, 1]
Ipwse    <- IpwEst[2, 2]
IpwCI    <- c(Ipwmean - 1.96*Ipwse, Ipwmean + 1.96*Ipwse)

#---DR
n      <- nrow(lalonde)
DRadj  <- 1/n^2 * (sum(sqrt((1 - PS1)/PS1)*(Linear1$fitted - mean(Y1))) + 
	sum(sqrt(PS0/(1 - PS0))*(Linear0$fitted - mean(Y0))))^2
DRse   <- sqrt(Ipwse^2 - DRadj)
DRmean <- DR1 - DR0
DRCI   <- c(DRmean - 1.96*DRse, DRmean + 1.96*DRse)


#---プロット
x    <- 1:4
mean <- c(Ttestmean, Psmean, Ipwmean, DRmean)
lcl  <- c(TtestCI[1], PsRestCI[1], IpwCI[1], DRCI[1])
ucl  <- c(TtestCI[2], PsRestCI[2], IpwCI[2], DRCI[2])

plot(x, mean, xaxt="n", ylim = c(0, 3500))
arrows(x, ucl, x, lcl, length=.05, angle=90, code=3) 
axis(side=1, at=1:4, labels=c("Ttest", "PS", "IPW", "DR"), cex=0.7)

f:id:isseing333:20110511231644j:image

このように、今回のデータでは調整してもしなくてもあまり結果は変わっていません(DRの計算は少し自信がありませんが、、)。


しかし、sparseなデータ(疎なデータ)では効果が逆転したり過剰に大きくなったりする事が多いです。

結果を良く吟味して利用すると良いかもしれません。


追記

多変量回帰と傾向スコアでの調整を比較した論文がありました。

http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B7P72-52W3V5B-2&_user=10&_coverDate=05%2F16%2F2011&_rdoc=1&_fmt=high&_orig=gateway&_origin=gateway&_sort=d&_docanchor=&view=c&_acct=C000050221&_version=1&_urlVersion=0&_userid=10&md5=5d4425d80a3e20b79c70783597dfd51f&searchtype=a

本当に傾向スコアって有用なのか?というタイトルで興味深いです。



追記2

傾向スコアの原論文はこれです。

http://faculty.smu.edu/millimet/classes/eco7377/papers/rosenbaum%20rubin%2083a.pdf

2011-05-02

小技メモ(Rで日本語の扱い)

Rは結構柔軟な言語になっていて、日本語もそのまま使えたりします。

例えばこんな感じでオブジェクト名を日本語にすることもできます。


> あ <- 1:10
> あ
 [1]  1  2  3  4  5  6  7  8  9 10

「あ」というオブジェクトに1〜10のベクトルを代入しています。

こんな感じでデータフレームの変数名も日本語に出来るのですが、調子に乗ってると痛い目に遭います。



実際に自分が痛い目を見たのですが、パッケージを作るときは日本語は受け付けません



日本語=不正なマルチバイト文字


です。



変数はアルファベットを使って指定し直せばいいのですが、例えば次のような時はどうしても日本語を使う必要があります。


  1. データに不正なマルチバイト文字が使われていて、パッケージの中で置換したい
  2. データ中の不正なマルチバイト文字で条件分岐させたい

これで私は詰まっていたのですが、twitterで@phosphor_mさんと@kos59125さんに相談に乗って頂いたおかげでなんとかなりました。

raw関数を使えば良くて、使用例は次のような感じです。


> as.numeric(charToRaw("あ"))
[1] 130 160
> rawToChar(as.raw(c(130, 160)))
[1] "あ"

このように、「あ」という文字列はRの中ではc(130, 160)というベクトルで認識されているようです。



もちろん次のような文字列もいけます。

> as.numeric(charToRaw("(´・ω・`) ショボーン"))
 [1]  40 129  76 129  69 131 214 129  69  96  41  32 131  86 131 135 131
[18] 123 129  91 131 147
> rawToChar(as.raw(c(40, 129,  76, 129,  69, 131, 214, 129,  69,  96, 
                    41,  32, 131,  86, 131, 135, 131, 123, 129,  91, 131, 147)))
[1] "(´・ω・`) ショボーン"
> 

これを利用することでパッケージ内で日本語を扱うことができました!

@phosphor_mさん、@kos59125さん、情報を下さってありあとうございます。



追記

この文字変換はSJIS環境でのコードなようです。

UTF-8環境では別の文字コードになるようでして、それはどこかにオプションを入れるのかな?

2011-05-01

小技メモ(macで日本語データを読み込む)

SJISはいつもイラっとさせられますが、mac版Rでデータが文字化けして困っている方も多いのでは。


このオプションを実行するだけで文字化けとはおさらば!(のはず笑)

options(encoding="SJIS")

だめだったらご一報下されば幸いです。