1. はじめに

みなさま、こんにちは。
冬ですね。街のイルミネーションや陽気な音楽が淋しくさせる季節になりましたね。切ないですね。
R Advent Calendar 2014 の7日目を担当させていただきます、aich_08_です。
どうぞよろしくお願い致します。

さて、突然ですが、

「●●推定量が漸近的に▲▲分布に従うことを数式を用いないで示してくれないかな」

と言われたことはありませんか？

データ解析をした結果を、社内であれば経営層、コンサルタントであればクライアントに報告する機会がありますが、
伝える相手は様々なバックグラウンドをもった方々で、そんな方にも解析に用いた手法を相手が理解できる言葉や方法で理解していただく必要がありますよね。

そんな時の一つの方法として、体感して理解していただくためのアプリケーションの例として、二項ロジットモデルにおけるパラメータの最尤推定値のヒストグラムと理論上の分布を比較するものを作ってみましたので、紹介します。

2. 最尤推定量の漸近的性質

数値シミュレーションに入る前に、簡単に最尤推定量の漸近的性質のおさらいです。
最尤推定量には、標本サイズが十分に大きくなった場合に真の値に確率収束するという一致性：$ \hat{\beta}^{MLE} \xrightarrow{p} \beta_0$や、
ある条件の下では、正規分布に分布収束するという漸近正規性：$\sqrt{n}\left(\hat{\beta}^{MLE} - \beta_0\right) \leadsto N(0,I(\beta_0)^{-1})$などがありますね。またこれらの性質を利用した様々な検定（尤度比検定、スコア検定、Wald検定）がありますね。漸近正規性は、対数尤度関数の$\beta_0$周りでのテイラー展開や、大数の法則、中心極限定理等を用いて証明されていて、なかなか難しいですよね。

そこで以下では上記に示す理論上の分布と数値シミュレーションにより得られる推定値の分布を比較する方法をまとめていきます。

3. 数値シミュレーション

今回は、次に示す二項ロジットモデルを考えます。

3.1. 設定

モデル : 説明変数$x_i$が与えられた下での、$0$または$1$の値を取る二値応答変数$y_i$は確率$p_i$で1をとる。$p_i$と$x_i$の間には$logit(p_i)=log\left(\frac{p_i}{1-p_i}\right) = \beta_0 x_i$の関係があるとする。
試行回数： $Nsim = 1000$
標本サイズ： $n = 250,500,750,1000$
説明変数 $x$ ：一様乱数($-10\leq x \leq 10 $)により生成
真のパラメータ値： $\beta_0 = 2.5$

3.2. シミュレーションの実行

以下には、Rで数値シミュレーションを実行するコードを示します。

準備

# 乱数のシード設定
set.seed(12345)

# 作業ディレクトリの設定
setwd("C:/hogehoge/R Advent Calendar")

# 真のパラメータ値を設定
beta1.o <- 2.5

試行1回を行うための関数を定義 : sim(n), where n : sample size.


# 試行1回分の関数
sim <- function(n){

  # 一様乱数を発生
  x1vec <- runif(n = n,min = -10,max = 10)
  zvec <- beta1.o*x1vec

  #p_iを計算
  pvec <- 1/(1+exp(-zvec))

  # 二項分布B(1,p_i)から応答変数y_iを求める
  yvec <- sapply(pvec,function(p) rbinom(1,1,p))

  sample.data <- data.frame(x=x1vec,y=yvec)

  # パラメータの推定
  res <- glm(y~x-1,family=binomial(link=logit),data=sample.data)

  # 最尤推定値
  est <- res$coefficients
  # 最適化が収束したか否か（0:収束していない,1:収束した）
  res.converged <- res$converged

  # 返り値の指定
  return(list(n=n,true=beta1.o,est=est,res.converged=res.converged))
}

シミュレーションの実行

# シミュレーションの設定
n.vec <- c(250,500,750,1000) #サンプルサイズ
Nsim <- 1000 # 試行回数

# シミュレーションの実行
for(i in seq(along.with=n.vec)){
  n <- n.vec[i]
  res.sim <- replicate(n = Nsim,sim(n))
  res.sim.mat <- matrix(unlist(res.sim),byrow=T,Nsim)
  if(i==1){
    res.all <- res.sim.mat
    colnames(res.all) <- c("n","true.b1","est.b1","converged")
  }else{
    res.all <- rbind(res.all,res.sim.mat)
  }
}

res.all.converged <- res.all[res.all[,"converged"]==1,]
write.csv(res.all.converged,"simulation.csv")

3.3. 可視化

ここからは、得られた推定値の分布をサンプルサイズ別に可視化し、理論上の分布の形状と比較するためのアプリケーションをつくります。
可視化には、shinyを用いますが、Rstudioのshinyチュートリアルページを参考にコーディングしました。http://shiny.rstudio.com/tutorial/

目指す形は、
* サンプルサイズ別の推定値のヒストグラム
* 真値の位置
* 理論上の漸近分布
が示されることです。

が、shinyの詳細をまだ理解できていないので、可視化までの手順とui.Rとserver.Rの中身だけ紹介します。

手順

Step.1 ：作業ディレクトリ直下に"App"というフォルダを作成する
Step.2 ： Appフォルダに以下に示すui.Rファイルおよびserver.Rファイルを格納する
- ui.Rファイル

library(shiny)

# Define UI for application that draws a histogram
shinyUI(fluidPage(

  # Application title
  titlePanel("MLE"),

  # Sidebar with a slider input for the sample size
  sidebarLayout(
    sidebarPanel(
      sliderInput("size",
                  "Sample size:",
                  min = 250,
                  max = 1000,
                  value = 250,
                  step = 250
                  )
    ),

    # Show a plot of the generated distribution
    mainPanel(
      plotOutput("distPlot")
    )
  )
))

server.Rファイル ※ 途中登場するフィッシャー情報量$I(\beta_0)$は別で計算しました。

library(shiny)
data <- read.csv("data/simulation.csv",header=TRUE)

# Define server logic required to draw a histogram
shinyServer(function(input, output) {

  output$distPlot <- renderPlot({
    x    <- subset(data,n==input$size,select=est.b1,drop = TRUE)  # Old Faithful Geyser data
    x.all <- subset(data,select=est.b1,drop=TRUE)
    bins <- 1 + log2(length(x))

    # 真値を与えた場合のフィッシャー情報量
    I0 <- 0.0105276

    # 理論上の分散
    t.variance <- 1/(input$size*I0)

    # 理論上の分布
    t.dist <- function(x){
      dnorm(x,mean = data$true.b1[1],sd = sqrt(t.variance))
    }


    # 推定値のヒストグラム
    hist(x, breaks = bins, col = 'darkgray', border = 'white',
         xlim=range(x.all),ylim=c(0,2),freq=FALSE,
         main = substitute(paste(label,hat(beta)^MLE),list(label="Histogram of ")),
         xlab = expression(hat(beta)^MLE))

    # 真値
    abline(v=data$true.b1[1],col="red")

    # 理論上の漸近分布
    curve(expr = t.dist,add = TRUE,col="blue")
    legend("topright",legend = c("true value","theoretical distribution"),lty=1,col=c("red","blue"))
  })
})

Step.3 ： Appフォルダ直下に"data"フォルダを作成し、数値シミュレーションの結果得られた"simulation.csv"を格納する
Step.4 ：アプリケーションを実行する

library(shiny)
runApp("App")

4. まとめ

二項ロジットモデルにおけるパラメータの最尤推定値を数値シミュレーションにより求め、その分布を確認しました。可視化には、当初ggvisを用いようと悪戦苦闘していたが、motion chartにはヒストグラムはない（？）。
見つからなかったのでshinyを用いて手動でサンプルサイズを変更するアプリケーションを作りました。
自動でサンプルサイズを増加させるアプリケーションを今後は作ってみたいと思います。