Se der certo teste esse exemplo.
library(rstan)
treatment_level <- c("0", "5", "10", "15")
repetition_number <- 20
response_mean <- c(100, 94, 88, 82)
treatment <- gl(n = length(treatment_level), k = repetition_number, labels = treatment_level)
str_error <- 2
response_temp <- sapply(response_mean, function(x){
response <- rnorm(n = repetition_number, mean = x, sd = str_error)
})
response <- c(response_temp)
plot(response ~ treatment)
summary(lm(response ~ treatment -1))
mean(response)
sd(response)
stan_data <- list(N = length(response), I = length(treatment_level), Treatment=as.numeric(treatment), Response = response)
# Bayesian method for estimate mean and standard error by treatment
stan_modelcode <- "
data { // data setup
int<lower=0> N; // sample size
int<lower=1> I; // number of treatments
real Response[N]; // Response
int<lower=1, upper=I> Treatment[N]; // Treatment
}
parameters {
real mu[I];
real<lower=0, upper=100> sigma[I];
}
model {
//Priors
mu ~ normal(0, 100);
sigma ~ uniform(0, 100);
//Likelihood
for(i in 1:N) {
Response[i] ~ normal(mu[Treatment[i]], sigma[Treatment[i]]);
}
}
"
fit <- stan(model_code = stan_modelcode, data = stan_data, iter = 1000, chains = 3,
verbose = TRUE)
plot(fit)