【问题标题】:Logistic regression in R, StanR,Stan中的逻辑回归
【发布时间】:2021-11-13 09:44:20
【问题描述】:

我想提取 Stan 拟合的预测值(在生成量块中)并将它们与实际观察结果进行比较,但我找不到简单的解决方案。下面是一个简单的逻辑回归模型是如何做到的:

library(rstan)
library(tidyverse)
library(boot)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

T <- 40
set.seed(123)
x <- sort(runif(T, 0, 10))
alpha <- 1
beta <- 0.2
logit_p <- alpha + beta * x
p <- inv.logit(logit_p)
y <- rbinom(T, 1, p)

model_code <- "
data {
  int<lower=0> N;
  vector[N] x;
  
  int<lower=0,upper=1> y[N];
  
}
parameters {
  real alpha;
  real beta;
}
model {
  y ~ bernoulli_logit(alpha + beta * x);
}
generated quantities {
  vector[N] z;
  for (n in 1:N)
    z[n] = bernoulli_logit_rng(alpha + beta * x[n]);
    
}"

model_data <- list(
  N = T, 
  x = x,
  y = y
)


stan_run <- stan(
  data = model_data,
  model_code = model_code
)


posterior <- rstan::extract(stan_run)

df <- as.data.frame(posterior$z)
df <-df %>% summarise(across(everything(.), 
                             ~ ifelse(length(.[which(. == 1)]) > length(.[which(. == 0)]), 1, 0)))

我不知道我的方法是否正确。有谁知道任何直接的方法吗?

【问题讨论】:

  • 我在分配 model_data 对象时收到错误“对象 x 未找到”。 xy 是什么?
  • @scrameri 我只是想展示我在拟合后提取后验的方法,因此 x 和 y 是什么无关紧要,您可以模拟您自己选择的任何 x 和 y 并将它们提供给模型。
  • 好的,但即使在模拟 x 和 y 之后我也无法拟合模型(语法错误:分配中的尺寸不匹配,第 18 行)。您的问题似乎更多是关于如何从 data.frame 中提取和总结某些东西,所以也许最好提供一个示例 posterior 对象,或者提供一个使用 rstan 和模拟 x 和 y 的可工作的、可重现的示例:stackoverflow.com/questions/5963269/…跨度>
  • @scrameri 我假设任何知道答案的人甚至根本不需要运行代码(我给出了多余的代码),因为这个问题非常笼统,与模型或数据(可以是任何数据或模型)。但是现在你问我编辑了代码,现在它应该可以工作了。
  • 是的,但是需要知道posterior 的结构是什么样的,这就是为什么现在有这个可重现的例子很好。不过,dput posterior 对象的一小部分就足够了。

标签: r stan


【解决方案1】:

apply 函数可以成为你的朋友:

运行示例

library(rstan)
#> Loading required package: StanHeaders
#> Loading required package: ggplot2
#> rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
#> For execution on a local, multicore CPU with excess RAM we recommend calling
#> options(mc.cores = parallel::detectCores()).
#> To avoid recompilation of unchanged Stan programs, we recommend calling
#> rstan_options(auto_write = TRUE)
library(tidyverse)
library(boot)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

T <- 40
set.seed(123)
x <- sort(runif(T, 0, 10))
alpha <- 1
beta <- 0.2
logit_p <- alpha + beta * x
p <- inv.logit(logit_p)
y <- rbinom(T, 1, p)

model_code <- "
data {
  int<lower=0> N;
  vector[N] x;
  
  int<lower=0,upper=1> y[N];
  
}
parameters {
  real alpha;
  real beta;
}
model {
  y ~ bernoulli_logit(alpha + beta * x);
}
generated quantities {
  vector[N] z;
  for (n in 1:N)
    z[n] = bernoulli_logit_rng(alpha + beta * x[n]);
    
}"

model_data <- list(
  N = T, 
  x = x,
  y = y
)

# run model
stan_run <- stan(
  data = model_data,
  model_code = model_code
)
#> Trying to compile a simple C file
#> Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
#> clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
#> In file included from <built-in>:1:
#> In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
#> In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
#> /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
#> namespace Eigen {
#> ^
#> /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
#> namespace Eigen {
#>                ^
#>                ;
#> In file included from <built-in>:1:
#> In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
#> /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#> #include <complex>
#>          ^~~~~~~~~
#> 3 errors generated.
#> make: *** [foo.o] Error 1

使用apply

在逻辑回归的情况下,您只需应用median 即可获得模态(最频繁)值。

(df.pred <- apply(rstan::extract(stan_run)$z, 2, median)

#>  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> [39] 1 1

reprex package (v2.0.1) 于 2021-09-19 创建

【讨论】:

  • 谢谢伙计,但这个解决方案与我的并没有太大不同(尽管它更好)。我有点想在 Stan 中找到一个可用的单线解决方案/功能。在 JAGS 中,您可以通过“model_fit$BUGSoutput$mean$z”来完成,因为 JAGS 为您完成,您只需访问它即可。我希望 Stan 也能这么简单。
  • 好吧,我认为这更多是为了简化您的 summarise 代码。我只是改变了我的解决方案来制作一个单线。我不知道有任何rstan 函数可以为您提取和总结它,但您可以尝试检查rstanarmrstantools 包。
猜你喜欢
  • 2018-01-26
  • 2018-02-12
  • 2014-06-20
  • 1970-01-01
  • 2014-06-26
  • 2019-08-20
  • 2021-05-02
  • 2017-07-30
  • 1970-01-01
相关资源
最近更新 更多