# Para instalar, descomente a linha abaixo
#install.packages("rstan", repos = c('https://stan-dev.r-universe.dev', getOption("repos")))Aula de Stan
Uma aula sobre Stan em R
Stan
Stan é uma linguagem de programação probabilística voltada para inferência estatística.
Seu nome é uma homenagem a Stanislaw Ulam, pioneiro nos métodos de Monte Carlo, conhecido também por trabalhar no projeto Manhattan, responsável pela morte de 130.000 a 250.000 pessoas, e por supostamente ter desenvolvido o primeiro projeto de bombas de hidrogênio.
O Stan foi desenvolvido por uma equipe de 52 pessoas e continua sendo atualizado constantemente.
Ele viabiliza modelos Bayesianos sofisticados, desde regressão linear simples até modelos multi-level e séries temporais.
Sua versatilidade se deve tanto a possibilidade de programar através de interfaces de outras linguagens como R, Pyhon e Julia, quanto a implementação dos algoritmos Monte Carlo Hamiltoniano (HMC) e No-U-Turn sampler (NUTS), além de permitir estimação de máxima verossimilhança (usualmente ignorada) e aproximações da posteriori (inferência variacional).
Vantagem do HMC em relação aos algoritmos mais simples como Metropolis-Hastings e amostrador de Gibbs se deve a acelerada convergência do processo em casos de altas dimensões pois gera passos longos com altas taxas de aceitação, reduzindo a correlação dos pontos.
Stan tem como concorrente o software Pymc, com este sendo exclusivo para Python.
Aprendendo Stan
Stan apresenta 3 documentos essenciais.
Stan Reference Manual: Manual de referência do usuário, válido para todas as interfaces da linguagem. A primeira parte consiste na descrição completa da linguagem, a segunda detalha os algoritmos e as ferramentas de inferência a posteriori, enquanto a terceira providencia informações axiliares do uso do Stan.
Stan User’s Guide: Guia do usuário. Contém modelos de exemplo e técnicas de programação.
Parte 1: Exibe códigos e discussões para várias classes de modelos.
Parte 2: Discute técnicas de programação
Parte 3: Apresenta algoritmos de calibração e verificação de modelos
Apêndice: Introdução do compilador, guia de estilo e dicas para usuários de outras ferramentas como BUGS (Bayesian Inference Using Gibbs Sampling) e JAGS (Just Another Gibbs Sampler)
Stan Functions Reference: Referência das funções implementadas.
Contudo, alguns materiais voltados para o ensino da linguagem estão disponíveis no próprio site, com destaque para o tutorial Getting Started with Bayesian Statistics using Stan and Python de Bob Carpenter (feito em Python). Além disso, vários estudos de casos também estão disponíveis como exemplos.
Os Pacotes
Como mencionado, Stan não é apenas um pacote, mas sim uma linguagem de programação por si só. Devido ao isso, para programar em Stan, é necessário o uso de uma interface. Como estamos usando R, duas interfaces estão disponíveis via pacotes.
CmdStanR: Interface mais leve ligada a versão mais recente do Stan. Simplesmente acessa o Stan e exibe as saídas. Originalmente mais limitada, mas o uso do pacote bridgestan ameniza tal disparidade. Indisponivel no CRAN
RStan: Interface original. Conecta-se diretamente com o código fonte do Stan, permitindo uma grande customização, mas atualizações são limitadas devido as restrições impostas pelo CRAN. Será o foco deste documento.
Além das interfaces, alguns pacotes auxiliares oficiais também estão disponiveis:
bayesplot: Baseado em ggplot2, apresenta diversos gráficos de amostragem da posteriori, diagnósticos visuais do MCMC e da distribuição preditiva a posteriori (ou priori)
BridgeStan: Pacote para acessar os métodos do Stan de forma eficiente e econômica
brms: Interface que adiciona modelos multilevel multivariado (não-)lineares com sintaxe similar ao pacote lme4, visando oferecer uma interface simples e familiar
loo: Implementa aproximação de leave-one-out cross-calibration para os modelos do Stan
posterior: Coleção de ferramentas auxiliares para análise da posteriori (ou priori) contribuindo para conversão de formato das saídas, operações comuns, medidas resumo e diagnóstico.
projpred: Implementa método para seleção de variáveis preditivas
rstanram: Pacote que emula outros pacotes de estatística do R, mas usando o Stan como base
rstantools: Ferramenta para a criação de pacotes
shinystan: Apresenta gráficos de diagnósticos interativos via shiny
Instalação
Para informações de como instalar, consulte https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started
Verificação da instalação:
# Para testar, descomente a linha abaixo
# example(stan_model, package = "rstan", run.dontrun = TRUE)Modelo linear
Primeiro, vamos chamar o pacote básico.
library("rstan") Loading required package: StanHeaders
rstan version 2.36.0.9000 (Stan version 2.37.0)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)
For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
change `threads_per_chain` option:
rstan_options(threads_per_chain = 1)# Se estiver numa máquina local, descomente as linhas abaixo para melhorar o desempenho
# options(mc.cores = parallel::detectCores())
# rstan_options(auto_write = TRUE)Vamos rodar um exemplo introdutório. Neste exemplo, estamos interesados em fazer a seguinte regressão linear nos dados cars
\[ \begin{aligned}y \mid \beta, \sigma^2 &\sim \mathcal{N}\bigl(X\beta,\, \sigma^2 I_N\bigr), \\[4pt]\beta \mid \sigma^2 &\sim \mathcal{N}\bigl(0,\, \sigma^2 \times 100\, I_J\bigr), \\[4pt]\sigma^2 &\sim \operatorname{InvGamma}(3,\,100),\end{aligned} \]
onde \(j=2\) é o número de parâmetros e \(N=50\) o número de indivíduos.
x = cars$dist # variável resposta
y = cars$speed # variável explicativa
n = length(x) # n=50
X = cbind(1,x) # Matrix de planejamento
p = ncol(X) # p=2
summary(X) V1 x
Min. :1 Min. : 2.00
1st Qu.:1 1st Qu.: 26.00
Median :1 Median : 36.00
Mean :1 Mean : 42.98
3rd Qu.:1 3rd Qu.: 56.00
Max. :1 Max. :120.00 Para definir o modelo, precisamos especificá-lo atráves da linguagem Stan e passar o arquivo para o modelo. Isso pode ser feito de duas maneiras, ou definindo um arquivo externo ou criando um arquivo de texto no R. De forma geral, o código é divido em três ou mais partes. Comumentemente, a inicialização das variáveis com as restrições, os parâmetros e o modelo em si.
rs_code <- '
data {
int<lower=1> N;
int<lower=1> J;
matrix[N,J] x;
vector[N] y;
}
parameters {
vector[J] beta;
real<lower=0> sigma2;
}
model {
sigma2 ~ inv_gamma(3, 100);
beta ~ normal(0, sqrt(sigma2*100));
y ~ normal(x * beta, sqrt(sigma2));
}'Como os ajustes são baseados em monte carlo, vamos definir o númeo de iterações, o tamanho do pulo da amostra a ser descartado e a amostra inicial a ser descartada. Com isso, rodamos o modelo.
burnin <- 2000
thin <- 3
N=(2000+burnin)*thin
# Conjunto de dados
stan_data <- list(N = n, J = p, y = y, x = X)
stan_mod <- stan(model_code = rs_code, data = stan_data,
chains = 4, iter = N, warmup = burnin, thin = thin)Trying to compile a simple C fileRunning /usr/local/lib/R/bin/R CMD SHLIB foo.c
using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
gcc -I"/usr/local/lib/R/include" -DNDEBUG -I"/usr/local/lib/R/site-library/Rcpp/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/unsupported" -I"/usr/local/lib/R/site-library/BH/include" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/src/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppParallel/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fpic -g -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g -c foo.c -o foo.o
In file included from <command-line>:
/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:3:10: fatal error: stdexcept: No such file or directory
3 | #include <stdexcept>
| ^~~~~~~~~~~
compilation terminated.
make: *** [/usr/local/lib/R/etc/Makeconf:195: foo.o] Error 1
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 1e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 12000 [ 0%] (Warmup)
Chain 1: Iteration: 1200 / 12000 [ 10%] (Warmup)
Chain 1: Iteration: 2001 / 12000 [ 16%] (Sampling)
Chain 1: Iteration: 3200 / 12000 [ 26%] (Sampling)
Chain 1: Iteration: 4400 / 12000 [ 36%] (Sampling)
Chain 1: Iteration: 5600 / 12000 [ 46%] (Sampling)
Chain 1: Iteration: 6800 / 12000 [ 56%] (Sampling)
Chain 1: Iteration: 8000 / 12000 [ 66%] (Sampling)
Chain 1: Iteration: 9200 / 12000 [ 76%] (Sampling)
Chain 1: Iteration: 10400 / 12000 [ 86%] (Sampling)
Chain 1: Iteration: 11600 / 12000 [ 96%] (Sampling)
Chain 1: Iteration: 12000 / 12000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.047 seconds (Warm-up)
Chain 1: 0.212 seconds (Sampling)
Chain 1: 0.259 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 4e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 12000 [ 0%] (Warmup)
Chain 2: Iteration: 1200 / 12000 [ 10%] (Warmup)
Chain 2: Iteration: 2001 / 12000 [ 16%] (Sampling)
Chain 2: Iteration: 3200 / 12000 [ 26%] (Sampling)
Chain 2: Iteration: 4400 / 12000 [ 36%] (Sampling)
Chain 2: Iteration: 5600 / 12000 [ 46%] (Sampling)
Chain 2: Iteration: 6800 / 12000 [ 56%] (Sampling)
Chain 2: Iteration: 8000 / 12000 [ 66%] (Sampling)
Chain 2: Iteration: 9200 / 12000 [ 76%] (Sampling)
Chain 2: Iteration: 10400 / 12000 [ 86%] (Sampling)
Chain 2: Iteration: 11600 / 12000 [ 96%] (Sampling)
Chain 2: Iteration: 12000 / 12000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.038 seconds (Warm-up)
Chain 2: 0.174 seconds (Sampling)
Chain 2: 0.212 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 3e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 12000 [ 0%] (Warmup)
Chain 3: Iteration: 1200 / 12000 [ 10%] (Warmup)
Chain 3: Iteration: 2001 / 12000 [ 16%] (Sampling)
Chain 3: Iteration: 3200 / 12000 [ 26%] (Sampling)
Chain 3: Iteration: 4400 / 12000 [ 36%] (Sampling)
Chain 3: Iteration: 5600 / 12000 [ 46%] (Sampling)
Chain 3: Iteration: 6800 / 12000 [ 56%] (Sampling)
Chain 3: Iteration: 8000 / 12000 [ 66%] (Sampling)
Chain 3: Iteration: 9200 / 12000 [ 76%] (Sampling)
Chain 3: Iteration: 10400 / 12000 [ 86%] (Sampling)
Chain 3: Iteration: 11600 / 12000 [ 96%] (Sampling)
Chain 3: Iteration: 12000 / 12000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.045 seconds (Warm-up)
Chain 3: 0.222 seconds (Sampling)
Chain 3: 0.267 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 4e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 12000 [ 0%] (Warmup)
Chain 4: Iteration: 1200 / 12000 [ 10%] (Warmup)
Chain 4: Iteration: 2001 / 12000 [ 16%] (Sampling)
Chain 4: Iteration: 3200 / 12000 [ 26%] (Sampling)
Chain 4: Iteration: 4400 / 12000 [ 36%] (Sampling)
Chain 4: Iteration: 5600 / 12000 [ 46%] (Sampling)
Chain 4: Iteration: 6800 / 12000 [ 56%] (Sampling)
Chain 4: Iteration: 8000 / 12000 [ 66%] (Sampling)
Chain 4: Iteration: 9200 / 12000 [ 76%] (Sampling)
Chain 4: Iteration: 10400 / 12000 [ 86%] (Sampling)
Chain 4: Iteration: 11600 / 12000 [ 96%] (Sampling)
Chain 4: Iteration: 12000 / 12000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.047 seconds (Warm-up)
Chain 4: 0.203 seconds (Sampling)
Chain 4: 0.25 seconds (Total)
Chain 4: Com o modelo ajustado, sem o uso de pacote adicionais, podemos obter informações a respeito de estatísticas resumo dos parâmetros a posteriori. A estatística n_eff é uma estimativa da amostra efetiva para cada parâmetro.
print(stan_mod)Inference for Stan model: anon_model.
4 chains, each with iter=12000; warmup=2000; thin=3;
post-warmup draws per chain=3334, total post-warmup draws=13336.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
beta[1] 8.29 0.01 0.98 6.35 7.64 8.30 8.94 10.19 11269 1
beta[2] 0.17 0.00 0.02 0.13 0.15 0.17 0.18 0.20 11331 1
sigma2 12.57 0.02 2.49 8.55 10.83 12.27 13.99 18.24 12851 1
lp__ -106.48 0.01 1.26 -109.79 -107.05 -106.16 -105.56 -105.05 10901 1
Samples were drawn using NUTS(diag_e) at Tue Oct 7 15:06:46 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).plot(stan_mod)ci_level: 0.8 (80% intervals)outer_level: 0.95 (95% intervals)
Podemos usar a função traceplot para observar a série temporal gerada pelo monte carlo e verificar se houve convergência.
traceplot(stan_mod,inc_warmup = TRUE, nrow = 3)
Podemos usar a função pair para verificar se houve problema na convergência. No gráfico gerado, valores com aceitação abaixo da mediana da taxa de aceitação do MCMC estão abaixos da diagonal principal e os maiores estão acima.
pairs(stan_mod,pars = c("beta[1]","beta[2]", "sigma2", "lp__"), las = 1)Warning in par(usr): argument 1 does not name a graphical parameter
Warning in par(usr): argument 1 does not name a graphical parameter
Warning in par(usr): argument 1 does not name a graphical parameter
Warning in par(usr): argument 1 does not name a graphical parameter
Vamos agora refazer as análises de diagnóstico usando o pacote bayesplot.
library("bayesplot")This is bayesplot version 1.14.0- Online documentation and vignettes at mc-stan.org/bayesplot- bayesplot theme set to bayesplot::theme_default() * Does _not_ affect other ggplot2 plots * See ?bayesplot_theme_set for details on theme settinglibrary("ggplot2")
posterior <- as.matrix(stan_mod)
plot_title <- ggtitle("Posterior distributions",
"with medians and 80% intervals")
mcmc_areas(posterior,
pars = c("beta[1]"),
prob = 0.8) + plot_title
mcmc_areas(posterior,
pars = c("beta[2]"),
prob = 0.8) + plot_title
mcmc_areas(posterior,
pars = c("sigma2"),
prob = 0.8) + plot_title
posterior2 <- extract(stan_mod, inc_warmup = TRUE, permuted = FALSE)
color_scheme_set("mix-blue-pink")
p <- mcmc_trace(posterior2, pars = c("beta[1]","beta[2]", "sigma2"), n_warmup = 300,
facet_args = list(nrow = 2, labeller = label_parsed))
p + facet_text(size = 15)
color_scheme_set("darkgray")
mcmc_scatter(
as.matrix(stan_mod),
pars = c("beta[1]", "beta[2]"),
np = nuts_params(stan_mod),
np_style = scatter_style_np(div_color = "green", div_alpha = 0.8)
)
color_scheme_set("darkgray")
mcmc_scatter(
as.matrix(stan_mod),
pars = c("beta[1]", "sigma2"),
np = nuts_params(stan_mod),
np_style = scatter_style_np(div_color = "green", div_alpha = 0.8)
)
color_scheme_set("red")
np <- nuts_params(stan_mod)
mcmc_nuts_energy(np,bins=60) + ggtitle("NUTS Energy Diagnostic")
mcmc_acf(
as.matrix(stan_mod),
pars = c("beta[1]", "beta[2]","sigma2"),
)
mcmc_dens_overlay(stan_mod, pars=c("beta[1]", "beta[2]")) 
De forma geral, os gráficos indicam que todas as cadeias convergiram. Podemos então gerar a preditiva a posteriori.
#---------------------------------------------
# Extração da posterior
#---------------------------------------------
post <- rstan::extract(stan_mod)
beta <- post$beta # matriz [n_iter x 2]
#---------------------------------------------
# Geração das linhas preditas
#---------------------------------------------
x_seq <- seq(min(x), max(x), length.out = 100)
X_seq <- cbind(1, x_seq)
# Selecionar um subconjunto de amostras (para não sobrecarregar o gráfico)
n_draws <- 200
idx <- sample(1:nrow(beta), n_draws)
# Criar data frame com todas as linhas preditas
lines_df <- data.frame()
for (i in idx) {
y_line <- X_seq %*% beta[i, ]
lines_df <- rbind(
lines_df,
data.frame(x = x_seq, y = y_line, draw = as.factor(i))
)
}
#---------------------------------------------
# Gráfico: todas as linhas geradas
#---------------------------------------------
ggplot() +
geom_point(aes(x = x, y = y), color = "black", size = 2) +
geom_line(
data = lines_df,
aes(x = x, y = y, group = draw),
color = "skyblue",
alpha = 0.2
) +
geom_line(
data = aggregate(y ~ x, data = lines_df, mean),
aes(x = x, y = y),
color = "blue",
linewidth = 1
) +
labs(
title = "Retas Preditivas a Posteriori",
subtitle = "Cada reta é uma amostra a posteriori de β",
x = "Distância",
y = "Velocidade"
) +
theme_minimal()
Modelo hierárquico binomial
Na aula 7, foi visto o modelo modelo hierárquico binomial a seguir \[ \begin{aligned} y_i \mid \theta_i &\sim \operatorname{Bin}\bigl(n_i, \, \theta_i \bigr), \\[4pt] \theta_i \mid \alpha,\,\beta &\sim \operatorname{Beta}(\alpha,\,\beta), \\[4pt] (\alpha,\,\beta) &\sim h(\alpha,\,\beta) \propto (\alpha+\beta)^{-5/2}, \\[4pt] \end{aligned} \]
Vamos agora replicar a análise do exercício 4 da lista 3 usando o mesmo modelo. No exercício, tinhamos a relação entre o número de bicicletas e o número de veículos automotres observados durante o período de 1 hora em quatro diferentes quadras. Essas quadras foram classificadas pelo tipo de tráfego (baixo ou moderado) e a existência ou não de ciclovias.
Suponha que o número de bicicletas siga um modelo binomial com probabilidade \(\theta_j,\,j=1,2,3,4\) e tamanho de amostra dado pelo número total de veículos.
# ----------------------------
# 2.1 Preparar os dados da tabela
# ----------------------------
# Grupo 1: Baixo, Ciclovia = "Baixo Sim"
y1 <- c(16, 9, 10, 13, 19, 20, 18, 17, 35, 55)
n1 <- c(58, 90, 48, 57, 103, 57, 86, 112, 273, 64)
# Grupo 2: Baixo, Não
y2 <- c(12, 1, 2, 4, 9, 7, 9, 8)
n2 <- c(113, 18, 14, 44, 208, 67, 29, 154)
# Grupo 3: Moderado, Ciclovia
y3 <- c(8, 35, 31, 19, 38, 47, 44, 44, 29, 18)
n3 <- c(29, 415, 425, 42, 180, 675, 620, 437, 47, 462)
# Grupo 4: Moderado, Não
y4 <- c(10, 43, 5, 14, 58, 15, 0, 47, 51, 32)
n4 <- c(557,1258,499,601,1163,700,90,1093,1459,1086)
# concatenar
y <- as.integer(c(y1, y2, y3, y4))
n <- as.integer(c(n1, n2, n3, n4))
# construir gid (grupo id 1..4)
gid <- as.integer(c(
rep(1, length(y1)),
rep(2, length(y2)),
rep(3, length(y3)),
rep(4, length(y4))
))
I <- length(y)
J <- 4
stan_data <- list(
I = I,
J = J,
y = y,
n = n,
gid = gid
)Com isso, podemos definir o seguinte modelo no Stan.
rs_code <-"
data {
int<lower=1> I;
int<lower=1> J;
array[I] int<lower=0> y;
array[I] int<lower=0> n;
array[I] int<lower=1, upper=J> gid;
}
parameters {
vector<lower=0, upper=1>[J] theta;
real<lower=0> alpha;
real<lower=0> beta;
}
model {
for (i in 1:I)
y[i] ~ binomial(n[i], theta[gid[i]]);
theta ~ beta(alpha, beta);
target += -2.5 * log(alpha + beta);
}
generated quantities {
real mean_theta;
real sd_theta;
mean_theta = mean(theta); // média dos 4 θ_j
sd_theta = sqrt(mean(square(theta - mean_theta)));
}
"Esse código apresenta três particularidades. A primeira, o uso do for dentro do modelo é possível, simplificando a escrita de modelos hierárquicos. Segundo, o uso do target += -2.5 * log(alpha + beta). O target no código stan representa a log densidade da distribuição a posteriori. Portanto, quando temos uma priori imprópria, sem distribuições, adicionamos o termo relativo a ela diretamente na densidade. De forma geral, o código inteiro poderia ser descrito com essas somas em target, bastando apresentar a log densidade de cada uma das funções. A terceira é o uso do bloco generated quantities. Este bloco não costuma ser necessário. Ele está sendo utilizado pois as quantias expressas nele são de interesse. É uma forma de nomear algumas variáveis importantes para facilitar o uso delas posteriormente.
burnin <- 2000
thin <- 3
N=(1250+burnin)*thin
stan_mod <- stan(model_code = rs_code, data = stan_data,
chains = 4, iter = N, warmup = burnin, thin = thin)Trying to compile a simple C fileRunning /usr/local/lib/R/bin/R CMD SHLIB foo.c
using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
gcc -I"/usr/local/lib/R/include" -DNDEBUG -I"/usr/local/lib/R/site-library/Rcpp/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/unsupported" -I"/usr/local/lib/R/site-library/BH/include" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/src/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppParallel/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fpic -g -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g -c foo.c -o foo.o
In file included from <command-line>:
/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:3:10: fatal error: stdexcept: No such file or directory
3 | #include <stdexcept>
| ^~~~~~~~~~~
compilation terminated.
make: *** [/usr/local/lib/R/etc/Makeconf:195: foo.o] Error 1
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 1.1e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 9750 [ 0%] (Warmup)
Chain 1: Iteration: 975 / 9750 [ 10%] (Warmup)
Chain 1: Iteration: 1950 / 9750 [ 20%] (Warmup)
Chain 1: Iteration: 2001 / 9750 [ 20%] (Sampling)
Chain 1: Iteration: 2975 / 9750 [ 30%] (Sampling)
Chain 1: Iteration: 3950 / 9750 [ 40%] (Sampling)
Chain 1: Iteration: 4925 / 9750 [ 50%] (Sampling)
Chain 1: Iteration: 5900 / 9750 [ 60%] (Sampling)
Chain 1: Iteration: 6875 / 9750 [ 70%] (Sampling)
Chain 1: Iteration: 7850 / 9750 [ 80%] (Sampling)
Chain 1: Iteration: 8825 / 9750 [ 90%] (Sampling)
Chain 1: Iteration: 9750 / 9750 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.083 seconds (Warm-up)
Chain 1: 0.344 seconds (Sampling)
Chain 1: 0.427 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 5e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 9750 [ 0%] (Warmup)
Chain 2: Iteration: 975 / 9750 [ 10%] (Warmup)
Chain 2: Iteration: 1950 / 9750 [ 20%] (Warmup)
Chain 2: Iteration: 2001 / 9750 [ 20%] (Sampling)
Chain 2: Iteration: 2975 / 9750 [ 30%] (Sampling)
Chain 2: Iteration: 3950 / 9750 [ 40%] (Sampling)
Chain 2: Iteration: 4925 / 9750 [ 50%] (Sampling)
Chain 2: Iteration: 5900 / 9750 [ 60%] (Sampling)
Chain 2: Iteration: 6875 / 9750 [ 70%] (Sampling)
Chain 2: Iteration: 7850 / 9750 [ 80%] (Sampling)
Chain 2: Iteration: 8825 / 9750 [ 90%] (Sampling)
Chain 2: Iteration: 9750 / 9750 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.085 seconds (Warm-up)
Chain 2: 0.381 seconds (Sampling)
Chain 2: 0.466 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 5e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 9750 [ 0%] (Warmup)
Chain 3: Iteration: 975 / 9750 [ 10%] (Warmup)
Chain 3: Iteration: 1950 / 9750 [ 20%] (Warmup)
Chain 3: Iteration: 2001 / 9750 [ 20%] (Sampling)
Chain 3: Iteration: 2975 / 9750 [ 30%] (Sampling)
Chain 3: Iteration: 3950 / 9750 [ 40%] (Sampling)
Chain 3: Iteration: 4925 / 9750 [ 50%] (Sampling)
Chain 3: Iteration: 5900 / 9750 [ 60%] (Sampling)
Chain 3: Iteration: 6875 / 9750 [ 70%] (Sampling)
Chain 3: Iteration: 7850 / 9750 [ 80%] (Sampling)
Chain 3: Iteration: 8825 / 9750 [ 90%] (Sampling)
Chain 3: Iteration: 9750 / 9750 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.079 seconds (Warm-up)
Chain 3: 0.341 seconds (Sampling)
Chain 3: 0.42 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 5e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 9750 [ 0%] (Warmup)
Chain 4: Iteration: 975 / 9750 [ 10%] (Warmup)
Chain 4: Iteration: 1950 / 9750 [ 20%] (Warmup)
Chain 4: Iteration: 2001 / 9750 [ 20%] (Sampling)
Chain 4: Iteration: 2975 / 9750 [ 30%] (Sampling)
Chain 4: Iteration: 3950 / 9750 [ 40%] (Sampling)
Chain 4: Iteration: 4925 / 9750 [ 50%] (Sampling)
Chain 4: Iteration: 5900 / 9750 [ 60%] (Sampling)
Chain 4: Iteration: 6875 / 9750 [ 70%] (Sampling)
Chain 4: Iteration: 7850 / 9750 [ 80%] (Sampling)
Chain 4: Iteration: 8825 / 9750 [ 90%] (Sampling)
Chain 4: Iteration: 9750 / 9750 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.077 seconds (Warm-up)
Chain 4: 0.313 seconds (Sampling)
Chain 4: 0.39 seconds (Total)
Chain 4: Warning: There were 5 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Warning: Examine the pairs() plot to diagnose sampling problemsO código avisou que neste caso havia alguns pontos com problema na convergência. Vamos rodar os gráficos para analisar.
pairs(stan_mod,pars = c("alpha","beta", "lp__"), las = 1)Warning in par(usr): argument 1 does not name a graphical parameter
Warning in par(usr): argument 1 does not name a graphical parameter
Warning in par(usr): argument 1 does not name a graphical parameter
posterior <- as.matrix(stan_mod)
plot_title <- ggtitle("Posterior distributions",
"with medians and 80% intervals")
mcmc_areas(posterior,
pars = c("alpha"),
prob = 0.8) + plot_title
mcmc_areas(posterior,
pars = c("beta"),
prob = 0.8) + plot_title
mcmc_areas(posterior,
pars = c("theta[1]"),
prob = 0.8) + plot_title
posterior2 <- extract(stan_mod, inc_warmup = TRUE, permuted = FALSE)
color_scheme_set("mix-blue-pink")
p <- mcmc_trace(posterior2, pars = c("alpha","beta", "theta[1]"), n_warmup = 300,
facet_args = list(nrow = 2, labeller = label_parsed))
p + facet_text(size = 15)
color_scheme_set("darkgray")
mcmc_scatter(
as.matrix(stan_mod),
pars = c("alpha", "beta"),
np = nuts_params(stan_mod),
np_style = scatter_style_np(div_color = "green", div_alpha = 0.8)
)
color_scheme_set("darkgray")
mcmc_scatter(
as.matrix(stan_mod),
pars = c("alpha", "theta[1]"),
np = nuts_params(stan_mod),
np_style = scatter_style_np(div_color = "green", div_alpha = 0.8)
)
color_scheme_set("red")
np <- nuts_params(stan_mod)
mcmc_nuts_energy(np,bins=60) + ggtitle("NUTS Energy Diagnostic")
mcmc_acf(
as.matrix(stan_mod),
pars = c("alpha", "beta","theta[1]"),
)
Aparentemente, apesar de algumas divergências, o modelo convergiu. Com isso, podemos obter algumas estatísticas e o gráfico das distribuiçoes a posteriori da média geral e do desvio padrão
\(\mathbb{E}[\theta_j \mid \alpha,\,\beta]\text{ e } DP[\theta_j \mid \alpha,\,\beta].\)
print(stan_mod)Inference for Stan model: anon_model.
4 chains, each with iter=9750; warmup=2000; thin=3;
post-warmup draws per chain=2584, total post-warmup draws=10336.
mean se_mean sd 2.5% 25% 50% 75% 97.5%
theta[1] 0.22 0.00 0.01 0.20 0.21 0.22 0.23 0.25
theta[2] 0.08 0.00 0.01 0.06 0.07 0.08 0.09 0.10
theta[3] 0.09 0.00 0.01 0.08 0.09 0.09 0.10 0.10
theta[4] 0.03 0.00 0.00 0.03 0.03 0.03 0.03 0.04
alpha 1.33 0.01 1.04 0.19 0.59 1.04 1.78 4.09
beta 9.09 0.10 9.15 0.42 2.57 6.16 12.58 33.91
mean_theta 0.11 0.00 0.00 0.10 0.10 0.11 0.11 0.12
sd_theta 0.07 0.00 0.01 0.06 0.07 0.07 0.07 0.08
lp__ -2948.54 0.02 1.76 -2952.87 -2949.48 -2948.23 -2947.25 -2946.11
n_eff Rhat
theta[1] 10290 1
theta[2] 10356 1
theta[3] 10623 1
theta[4] 10037 1
alpha 9083 1
beta 9170 1
mean_theta 10761 1
sd_theta 10196 1
lp__ 8862 1
Samples were drawn using NUTS(diag_e) at Tue Oct 7 15:07:45 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).mcmc_areas(posterior,
pars = c("mean_theta"),
prob = 0.9) + plot_title
mcmc_areas(posterior,
pars = c("sd_theta"),
prob = 0.9) + plot_title
Comentários sobre os blocos do Stan
O arquivo .stan é organizado em blocos:
functions: define funções a serem utilizadas
data: declara os dados a serem utilizados
transformed data: define transfornações dos dados
parameters: declara os parâmetros a serem estimados
transformed parameters: define transformações dos parâmetros
model: especifica o modelo
generated quantities: define valores desejados na saída
Em geral, bastaria os arumentos data, parameters e model para fazer qualquer análise, mas os outros blocos permitem facilitar e agilizar o trabalho.
Modelo linear generalizado
Para fazer modelos da classe do MLG, bastaria apenas alterar o código colocando as distribuições desejadas como resposta. Contudo, o pacote rstanarm pode ser utilizado no lugar do código Stan, trazendo uma notação muito similar à utilizada em MLG. Vamos considerar o exemplo de McCullagh & Nelder (1989). Nele, estamos analisando o tempo até a coagulação sanguínea para dois grupos diferentes tendo como variável explicativa a concentração de plasma livre de protombina. O modelo usual frequentista seria o abaixo.
clotting <- data.frame(
u = c(5,10,15,20,30,40,60,80,100),
lot1 = c(118,58,42,35,27,25,21,19,18),
lot2 = c(69,35,26,21,18,16,13,12,12))
clotting2 <- with(clotting, data.frame(
log_plasma = rep(log(u), 2),
clot_time = c(lot1, lot2),
lot_id = factor(rep(c(1,2), each = length(u)))
))
summary(glm(clot_time ~ log_plasma * lot_id, data = clotting2, family = Gamma))
Call:
glm(formula = clot_time ~ log_plasma * lot_id, family = Gamma,
data = clotting2)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0165544 0.0008655 -19.127 1.97e-11 ***
log_plasma 0.0153431 0.0003872 39.626 8.85e-16 ***
lot_id2 -0.0073541 0.0016780 -4.383 0.000625 ***
log_plasma:lot_id2 0.0082561 0.0007353 11.228 2.18e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.002129707)
Null deviance: 7.708667 on 17 degrees of freedom
Residual deviance: 0.029401 on 14 degrees of freedom
AIC: 63.195
Number of Fisher Scoring iterations: 3Contudo, podemos facilmente replicar essa análise de forma bayesiana usando a função stan_glm.
library(rstanarm)Loading required package: RcppThis is rstanarm version 2.32.2- See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!- Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.- For execution on a local, multicore CPU with excess RAM we recommend calling options(mc.cores = parallel::detectCores())
Attaching package: 'rstanarm'The following object is masked from 'package:rstan':
loofit <- stan_glm(clot_time ~ log_plasma * lot_id, data = clotting2, family = Gamma,
prior_intercept = normal(0, 1, autoscale = TRUE),
prior = normal(0, 1, autoscale = TRUE))
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 4.1e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.461 seconds (Warm-up)
Chain 1: 0.295 seconds (Sampling)
Chain 1: 0.756 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 1e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.497 seconds (Warm-up)
Chain 2: 0.183 seconds (Sampling)
Chain 2: 0.68 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 1e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.456 seconds (Warm-up)
Chain 3: 0.097 seconds (Sampling)
Chain 3: 0.553 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 1e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.454 seconds (Warm-up)
Chain 4: 0.168 seconds (Sampling)
Chain 4: 0.622 seconds (Total)
Chain 4: print(fit,digits=3)stan_glm
family: Gamma [inverse]
formula: clot_time ~ log_plasma * lot_id
observations: 18
predictors: 4
------
Median MAD_SD
(Intercept) -0.016 0.007
log_plasma 0.015 0.003
lot_id2 -0.006 0.014
log_plasma:lot_id2 0.008 0.006
Auxiliary parameter(s):
Median MAD_SD
shape 7.631 2.739
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregE também, a mesma análise de convergência também pode ser realizada.
posterior <- as.matrix(fit)
plot_title <- ggtitle("Posterior distributions",
"with medians and 80% intervals")
mcmc_areas(posterior,
pars = c("(Intercept)"),
prob = 0.8) + plot_title
mcmc_areas(posterior,
pars = c("log_plasma"),
prob = 0.8) + plot_title
mcmc_areas(posterior,
pars = c("lot_id2"),
prob = 0.8) + plot_title
mcmc_areas(posterior,
pars = c("log_plasma:lot_id2"),
prob = 0.8) + plot_title
mcmc_areas(posterior,
pars = c("shape"),
prob = 0.8) + plot_title
color_scheme_set("darkgray")
mcmc_scatter(
as.matrix(fit),
pars = c("log_plasma","lot_id2"),
np = nuts_params(fit),
np_style = scatter_style_np(div_color = "green", div_alpha = 0.8)
)
color_scheme_set("red")
np <- nuts_params(fit)
mcmc_nuts_energy(np,bins=60) + ggtitle("NUTS Energy Diagnostic")
mcmc_acf(
as.matrix(fit),
pars = c("log_plasma","lot_id2"),
)
Uma alternative para o rstanarm é o pacote brms. Ambos os pacotes são simplificações do código para uma notação mais similar à usual dos pacotes estatísticos. Contudo, enquanto a primeira é mais rápida, ela é bastante limitada pois depende dos modelos já implementados, já o brms adiciona uma variadade maior de modelos, incluindo não lineares e aditivos. Abaixo replicamos o último modelo com a função brm.
library(brms)Loading 'brms' package (version 2.23.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: 'brms'The following objects are masked from 'package:rstanarm':
dirichlet, exponential, get_y, lasso, ngrpsThe following object is masked from 'package:bayesplot':
rhatThe following object is masked from 'package:rstan':
looThe following object is masked from 'package:stats':
arfit_brms <- brm(
clot_time ~ log_plasma * lot_id,
data = clotting2,
family = Gamma(link = "inverse"),
prior = c(
prior(normal(0, 1), class = "Intercept"),
prior(normal(0, 1), class = "b")
),
chains = 4, iter = 2000, warmup = 1000
)Compiling Stan program...Trying to compile a simple C fileRunning /usr/local/lib/R/bin/R CMD SHLIB foo.c
using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
gcc -I"/usr/local/lib/R/include" -DNDEBUG -I"/usr/local/lib/R/site-library/Rcpp/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/unsupported" -I"/usr/local/lib/R/site-library/BH/include" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/src/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppParallel/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fpic -g -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g -c foo.c -o foo.o
In file included from <command-line>:
/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:3:10: fatal error: stdexcept: No such file or directory
3 | #include <stdexcept>
| ^~~~~~~~~~~
compilation terminated.
make: *** [/usr/local/lib/R/etc/Makeconf:195: foo.o] Error 1Start sampling
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[9] is -1.06862, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[7] is -0.123557, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.12616, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[10] is -0.0168284, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[12] is -0.313332, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.11298, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.604841, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[14] is -1.25312, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[10] is -3.4122, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[1] is -27.675, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[9] is -0.0838248, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[12] is -0.0391665, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[2] is -0.254452, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.495122, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[6] is -0.0254338, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[8] is -0.00485711, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[12] is -0.456653, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: gamma_lpdf: Inverse scale parameter[6] is -0.104256, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 1:
Chain 1: Gradient evaluation took 1.5e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.759 seconds (Warm-up)
Chain 1: 0.113 seconds (Sampling)
Chain 1: 0.872 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -23.125, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[10] is -0.375688, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[5] is -0.260904, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.77424, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[10] is -4.07644, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[2] is -1.30262, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.10007, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[10] is -0.137297, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[2] is -1.49123, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[12] is -0.135085, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.89699, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.77452, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -2.37181, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -4.31924, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.6689, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.191961, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[8] is -0.0647898, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.444468, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.02738, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[8] is -1.56945, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.232961, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -17.2214, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.10355, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[4] is -0.50837, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[10] is -0.504332, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -13.2623, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -42.6027, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[4] is -0.0433483, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[6] is -0.055019, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.539601, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[11] is -0.00512004, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -5.88879, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.241846, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 2:
Chain 2: Gradient evaluation took 6e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.883 seconds (Warm-up)
Chain 2: 0.104 seconds (Sampling)
Chain 2: 0.987 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[6] is -0.0248506, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.717646, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[11] is -0.793113, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[2] is -1.34171, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[11] is -0.336961, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[10] is -1.06394, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[3] is -0.0257783, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.89696, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[12] is -1.11063, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -26.7572, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[10] is -0.555497, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.708408, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[12] is -0.0146853, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -3.94415, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -13.8846, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[5] is -0.0192127, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[13] is -0.0936162, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -2.98896, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.116842, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -3.6312, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -10.6328, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[10] is -7.8314, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.00474, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[10] is -0.0834285, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.33178, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.0957, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.59482, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -13.8761, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -11.5163, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[12] is -0.20574, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.820895, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -9.10549, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -21.3055, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[4] is -0.17643, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: gamma_lpdf: Inverse scale parameter[1] is -15.3485, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 3:
Chain 3: Gradient evaluation took 6e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.345 seconds (Warm-up)
Chain 3: 0.114 seconds (Sampling)
Chain 3: 0.459 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -2.9428, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[10] is -5.30925, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.353795, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -1.29255, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -14.9113, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[2] is -0.100189, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.233763, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[7] is -0.304611, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[12] is -2.83186, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.637665, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.915707, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.108764, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -4.00305, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -9.63312, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[2] is -0.343649, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[8] is -0.0609529, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -2.98937, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[10] is -0.682656, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[3] is -0.869155, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: gamma_lpdf: Inverse scale parameter[1] is -14.8845, but must be positive finite! (in 'anon_model', line 39, column 4 to column 49)
Chain 4:
Chain 4: Gradient evaluation took 6e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.345 seconds (Warm-up)
Chain 4: 0.119 seconds (Sampling)
Chain 4: 0.464 seconds (Total)
Chain 4: Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Warning: Examine the pairs() plot to diagnose sampling problemsOs resultados de todos os modelos concordam entre si.
print(fit_brms)Warning: There were 1 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup Family: gamma
Links: mu = inverse
Formula: clot_time ~ log_plasma * lot_id
Data: clotting2 (Number of observations: 18)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.02 0.00 -0.02 -0.01 1.00 2635 2429
log_plasma 0.02 0.00 0.01 0.02 1.00 2602 2248
lot_id2 -0.01 0.00 -0.01 -0.00 1.00 1383 1430
log_plasma:lot_id2 0.01 0.00 0.01 0.01 1.00 1491 1500
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 283.60 107.94 109.48 530.95 1.00 792 1148
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).plot(fit_brms)
Modelo aditivo generalizado
O pacote brms também permite o ajuste de modelos não linear. Como exemplo, vamos ajustar os dados de temperatura anual média do ar sobre o solo do planeta terra, disponível no pacote aida.
#install.packages('remotes')
#remotes::install_github('michael-franke/aida-package')
data_temperature <- aida::data_WorldTempggplot(data_temperature, aes(x = year, y = avg_temp)) +
geom_line(color = "steelblue", linewidth = 1) +
labs(
x = "Ano",
y = "Temperatura Média (°C)"
) +
theme_minimal(base_size = 14)
Podemos ajustar o modelo linear tradicional a seguir.
fit_temperature <- brm(formula = avg_temp ~ year,
data = data_temperature, refresh=0
)Compiling Stan program...Trying to compile a simple C fileRunning /usr/local/lib/R/bin/R CMD SHLIB foo.c
using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
gcc -I"/usr/local/lib/R/include" -DNDEBUG -I"/usr/local/lib/R/site-library/Rcpp/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/unsupported" -I"/usr/local/lib/R/site-library/BH/include" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/src/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppParallel/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fpic -g -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g -c foo.c -o foo.o
In file included from <command-line>:
/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:3:10: fatal error: stdexcept: No such file or directory
3 | #include <stdexcept>
| ^~~~~~~~~~~
compilation terminated.
make: *** [/usr/local/lib/R/etc/Makeconf:195: foo.o] Error 1Start samplinglibrary(dplyr)
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionpred_data <- data_temperature %>%
data.frame(year = .$year) %>%
mutate(pred = fitted(fit_temperature, newdata = ., probs = c(0.05, 0.95))[, "Estimate"],
lwr = fitted(fit_temperature, newdata = ., probs = c(0.05, 0.95))[, "Q5"],
upr = fitted(fit_temperature, newdata = ., probs = c(0.05, 0.95))[, "Q95"])
# Plot: observed data + posterior mean + 90% credible interval
ggplot(pred_data, aes(x = year, y = avg_temp)) +
geom_point(color = "gray40", alpha = 0.6) +
geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "skyblue", alpha = 0.3) +
geom_line(aes(y = pred), color = "firebrick", linewidth = 1.2) +
labs( x = "Ano",
y = "Temperatura Média"
) +
theme_minimal(base_size = 14)
Contudo, tal modelo não necessariamente é adequado visto que supõem crescimento linear para a temperatura. Devido a isso, um modelo alternativo, e talvez mais apropriado, seja baseado em spline para efeito do ano.
fit_temperature <- brm(formula = avg_temp ~ s(year),
data = data_temperature, refresh=0
)Compiling Stan program...Trying to compile a simple C fileRunning /usr/local/lib/R/bin/R CMD SHLIB foo.c
using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
gcc -I"/usr/local/lib/R/include" -DNDEBUG -I"/usr/local/lib/R/site-library/Rcpp/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/unsupported" -I"/usr/local/lib/R/site-library/BH/include" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/src/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppParallel/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fpic -g -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g -c foo.c -o foo.o
In file included from <command-line>:
/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:3:10: fatal error: stdexcept: No such file or directory
3 | #include <stdexcept>
| ^~~~~~~~~~~
compilation terminated.
make: *** [/usr/local/lib/R/etc/Makeconf:195: foo.o] Error 1Start samplingpred_data <- data_temperature %>%
data.frame(year = .$year) %>%
mutate(pred = fitted(fit_temperature, newdata = ., probs = c(0.05, 0.95))[, "Estimate"],
lwr = fitted(fit_temperature, newdata = ., probs = c(0.05, 0.95))[, "Q5"],
upr = fitted(fit_temperature, newdata = ., probs = c(0.05, 0.95))[, "Q95"])
# Plot: observed data + posterior mean + 90% credible interval
ggplot(pred_data, aes(x = year, y = avg_temp)) +
geom_point(color = "gray40", alpha = 0.6) +
geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "skyblue", alpha = 0.3) +
geom_line(aes(y = pred), color = "firebrick", linewidth = 1.2) +
labs( x = "Ano",
y = "Temperatura Média"
) +
theme_minimal(base_size = 14)
Análise de sobrevivência
Uma outra classe de modelos interessante é a de análise de sobrevivência. Nessa classe de modelos, usualmente estudamos o tempo até a ocorrência de eventos. Vamos analisar o tempo de recorrência de infecção em pacientes com doença renal.
library(survival)
Attaching package: 'survival'The following object is masked from 'package:brms':
kidneyhead(kidney, n = 5) id time status age sex disease frail
1 1 8 1 28 1 Other 2.3
2 1 16 1 28 1 Other 2.3
3 2 23 1 48 2 GN 1.9
4 2 13 0 48 2 GN 1.9
5 3 22 1 32 1 Other 1.2plot(survfit(Surv(time, status) ~ 1, data = kidney), xlab='dias', main = 'Kaplan-Meier')
Podemos então ajustar um modelo lognormal para o tempo até a recorrência. Inclusive, estaremos ajustando um modelo de efeito aleatório (também denominado de modelo de fragilidade no contexto de sobrevivência) e lkj é a distribuição a priori da matriz de correlações (Lewandowski, Kurowicka, and Joe, 2009).
fit1 <- brm(
formula = time | cens(status-1) ~ age * sex + disease + (1 + age|id),
data = kidney, family = lognormal(),
prior = c(set_prior("normal(0,5)", class = "b"),
set_prior("cauchy(0,2)", class = "sd"),
set_prior("lkj(2)", class = "cor")), warmup = 1000,
iter = 2000, chains = 4, control = list(adapt_delta = 0.95),
refresh=0)Compiling Stan program...Trying to compile a simple C fileRunning /usr/local/lib/R/bin/R CMD SHLIB foo.c
using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
gcc -I"/usr/local/lib/R/include" -DNDEBUG -I"/usr/local/lib/R/site-library/Rcpp/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppEigen/include/unsupported" -I"/usr/local/lib/R/site-library/BH/include" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/src/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/RcppParallel/include/" -I"/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fpic -g -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g -c foo.c -o foo.o
In file included from <command-line>:
/home/posmae/eduardojanotti/R/x86_64-pc-linux-gnu-library/4.4/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:3:10: fatal error: stdexcept: No such file or directory
3 | #include <stdexcept>
| ^~~~~~~~~~~
compilation terminated.
make: *** [/usr/local/lib/R/etc/Makeconf:195: foo.o] Error 1Start samplingWarning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Warning: Examine the pairs() plot to diagnose sampling problemsprint(fit1)Warning: There were 1 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup Family: lognormal
Links: mu = identity
Formula: time | cens(status - 1) ~ age * sex + disease + (1 + age | id)
Data: kidney (Number of observations: 76)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~id (Number of levels: 38)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.41 0.31 0.02 1.17 1.00 1777 2035
sd(age) 0.01 0.01 0.00 0.03 1.00 1264 2212
cor(Intercept,age) -0.12 0.46 -0.87 0.77 1.00 2670 2554
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.30 2.48 -2.53 7.14 1.00 1978 2591
age -0.00 0.06 -0.12 0.11 1.00 1872 2567
sex 1.01 1.40 -1.71 3.71 1.00 1943 2450
diseaseGN -0.26 0.66 -1.56 1.01 1.00 3224 3163
diseaseAN -0.31 0.64 -1.60 0.94 1.00 3108 3234
diseasePKD 0.51 0.89 -1.23 2.25 1.00 3502 3287
age:sex -0.00 0.03 -0.07 0.06 1.00 1842 2067
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.54 0.18 1.23 1.94 1.00 3235 3112
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).plot(fit1)
# Kaplan–Meier estimate
km <- survfit(Surv(time, status) ~ 1, data = kidney)
# Create a grid of times to evaluate survival
tseq <- seq(0, max(kidney$time), length.out = 200)
# Draw 100 posterior samples of linear predictor (no random effects)
lp <- posterior_linpred(fit1, transform = TRUE, re_formula = NA, ndraws = 100)
# Compute survival for each draw
# For lognormal, survival = 1 - Phi((log(t) - mu)/sigma)
sigma <- exp(posterior_samples(fit1, "sigma")[1:100, ])Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.mu <- lp[, 1] # first observation’s linear predictor
S <- sapply(1:100, function(i) {
1 - pnorm((log(tseq) - mu[i]) / sigma[i])
})
# Plot Kaplan–Meier
plot(km, conf.int = FALSE, lwd = 2, col = "black",
xlab = "Tempo", ylab = "Sobrevivência")
# Add posterior survival curves (semi-transparent blue)
matlines(tseq, S, col = rgb(0, 0, 1, 0.1), lty = 1)
Aparentemente, mesmo com a adição das covariáveis, o modelo não parece concordar com as curvas de Kaplan-Meier. Uma alternativa seria o uso do modelo de riscos proporcionais de Cox. Contudo, tal classe de modelos já está além do escorpo do presente material.
Previsão de séries temporais
Um outro pacote de destaque relacionado ao Stan é o pacote prophet. Ele foi criado pelos desenvolvedores da Meta (então facebook) para previsão de séries temporais baseado em Stan. Vamos fazer a previsão da série analisada anteriormente das temperaturas
library(prophet)Loading required package: rlangdata_temperature$ds <- as.Date(paste0(data_temperature$year, "-01-01"))
data_temperature$y <- data_temperature$avg_temp
df_for_model <- data_temperature[, c("ds", "y")]
df_for_model <- na.omit(df_for_model)
mod = prophet(df_for_model,weekly.seasonality=FALSE,daily.seasonality=FALSE)
future <- make_future_dataframe(mod,freq="year", periods = 100)
tail(future) ds
364 2114-01-01
365 2115-01-01
366 2116-01-01
367 2117-01-01
368 2118-01-01
369 2119-01-01forecast <- predict(mod, future)
tail(forecast[c('ds', 'yhat', 'yhat_lower', 'yhat_upper')]) ds yhat yhat_lower yhat_upper
364 2114-01-01 11.27551 10.74403 11.78980
365 2115-01-01 11.28803 10.75919 11.81828
366 2116-01-01 11.30593 10.74084 11.82309
367 2117-01-01 11.33980 10.78668 11.84780
368 2118-01-01 11.34693 10.82000 11.89496
369 2119-01-01 11.35945 10.82567 11.91653plot(mod, forecast)
prophet_plot_components(mod, forecast)