0

問題

小さなデータセット (N=100) があります。ポアソン回帰を実行する必要がありますが、一度に 1 つの観測を除外します (したがって、ローリングポアソン回帰)。

方程式にはいくつかの予測子がありますが、私が気にするのは 1 つです ( bxと呼びます)。私の考えは、bx が 100 個のモデル間でどの程度変化するかを見ることです。次に、これらの 100 ポイントの推定値を、Y 軸に効果の大きさ、X 軸にモデル番号を付けてプロットしたいと思います。

要約すると、次のものが必要です。

  1. JAGS でローリング ポアソン回帰を実行します (R2jags 経由)。

  2. 見積もりを取得したら、それらをプロットします。

JAGS の私のポアソン モデルは正常に動作していることに注意してください(以下は、私のモデル/データのサンプル おもちゃです)。ただし、「ローリング」バージョンを実装できていません。

自己完結型の例

# clear R
rm(list=ls())
cat("\014")

# load libraries
if (!require("pacman")) install.packages("pacman"); library(pacman) 
p_load(R2jags)


# Toy Data
N <- 100
x <- rnorm(n=N)  # standard Normal predictor
y <- rpois(n=N, lambda = 1)  # Poisson DV


# model
model <- function() {
    ## Likelihood
    for(i in 1:N){
            y[i] ~ dpois(lambda[i])
            log(lambda[i]) <- 
                    mu + # intercept
                    b.x*x[i]
            }
    ## Priors 
    mu ~  dnorm(0, 0.01) ## intercept
    b.x ~ dnorm(0, 0.01)
    }

# list elements
data.list <- list(N = N, y = y, x = x)

# run model
model.fit <- jags(
    data=data.list,
    inits=NULL,
    parameters.to.save = c("b.x"),
    n.chains = 1,
    n.iter = 20,
    n.burnin = 2, 
    model.file=model,
    progress.bar = "none")

Ok。それがモデルです。bxmodel.fitがあり、100回取得する必要がある係数です。現在のコードでは、完全なデータセットで一度だけ取得できます。ただし、df の最初の行を除外して 2 回目、次に df の 2 行目を除外して 3 回目、というように取得する必要があります。そして、これらすべてのbxをプロットします。

ここで、例として、最初の要素 ( bxの係数) が必要であることを知らせるために、単純なテーブルを作成します。

## I sourced this function below from    https://raw.githubusercontent.com/jkarreth/JKmisc/master/mcmctab.R

# Function to Create Table
mcmctab <- function(sims, ci = .8, digits = 2){

    require(coda) 

    if(class(sims) == "jags" | class(sims) == "rjags"){
            sims <- as.matrix(as.mcmc(sims))
    }
    if(class(sims) == "bugs"){
            sims <- sims$sims.matrix
    }  
    if(class(sims) == "mcmc"){
            sims <- as.matrix(sims)
    }    
    if(class(sims) == "mcmc.list"){
            sims <- as.matrix(sims)
    }      
    if(class(sims) == "stanfit"){
            stan_sims <- rstan::As.mcmc.list(sims)
            sims <- as.matrix(stan_sims)
    }      


    dat <- t(sims)

    mcmctab <- apply(dat, 1,
                     function(x) c(Mean = round(mean(x), digits = digits), # Posterior mean
                                   SD = round(sd(x), digits = 3), # Posterior SD
                                   Lower = as.numeric(
                                           round(quantile(x, probs = c((1 - ci) / 2)), 
                                                 digits = digits)), # Lower CI of posterior
                                   Upper = as.numeric(
                                           round(quantile(x, probs = c((1 + ci) / 2)), 
                                                 digits = digits)), # Upper CI of posterior
                                   Pr. = round(
                                           ifelse(mean(x) > 0, length(x[x > 0]) / length(x),
                                                  length(x[x < 0]) / length(x)), 
                                           digits = digits) # Probability of posterior >/< 0
                     ))
    return(t(mcmctab))
}


# this is the coefficient I need, but with different data frames.
mcmctab(model.fit)[1,1]

申し訳ありませんが、ここで試みられた解決策を提供することさえできません。よろしくお願いします。

4

2 に答える 2