I am new to the R stan language but I had an ode compartmental which I wrote in R stan. Now I want to integrate it and fit my actual and observed data. Any help will be appreciated. This is the stan file code and the stan data which I have worked on. The data are counts hence I started working using negative binomial due to overdispersion but the simulated data presents a poor fit to my data. I will also be glad if I can implement both the negative binomial and the ode integration and fit the observed and simulated data in R stan.
data {
int<lower=0> T; // Number of time steps
real times[T]; // Time steps
real rainfall[T]; // Rainfall data
real temperature[T]; // Temperature data
real bites[T]; // Biting rate data
real itn_cov[T]; // ITN coverage
real irs_cov[T]; // IRS coverage
real<lower=0> observed_ICH[T]; //observed uncomplicated malaria
real<lower=0> observed_Sev[T]; //observed severe malaria
}
parameters {
real<lower=0> SH;
real<lower=0> EH;
real<lower=0> ICH;
real<lower=0> IAH;
real<lower=0> Sev;
real<lower=0> TACTH;
real<lower=0> TASH;
real<lower=0> TFACTH;
real<lower=0> TFASH;
real<lower=0> AM;
real<lower=0> SM;
real<lower=0> EM;
real<lower=0> IM;
real<lower=0> m2;
real<lower=0> ah1;
real<lower=0> ah2;
real<lower=0> rs1;
real<lower=0> rs2;
real<lower=0> rho1;
real<lower=0> rho2;
real<lower=0> pc2;
real<lower=0> gamma;
real<lower=0> q;
real<lower=0> pt1;
real<lower=0> tau;
real<lower=0> pst;
real<lower=0> RDT;
real<lower=0> ps;
real<lower=0> ss;
real<lower=0> omega;
real<lower=0> pt2;
real<lower=0> trf1;
real<lower=0> trf2;
real<lower=0> Kc;
real<lower=0> betaMH;
real<lower=0> betaHM;
real<lower=0> itn_eff;
real<lower=0> itn_aversion;
real<lower=0> itn_use;
real<lower=0> irs_eff;
real<lower=0> irs_aversion;
// real<lower = 0.3, upper = 0.7> m2;
// real<lower = 0.3, upper = 0.5> ah1;
// real<lower = 0.05, upper = 0.15> ah2;
// real<lower = 0.02, upper = 0.06> rs1;
// real<lower = 0.005, upper = 0.015> rs2;
// real<lower = 365.25/4, upper = 365.25/2> rho1;
// real<lower = 365.25/7, upper = 365.25/5> rho2;
// real<lower = 0.1, upper = 0.2> pc2;
// real<lower = 365.25/25, upper = 365.25/17> gamma;
// real<lower = 365.25/220, upper = 365.25/170> q;
// real<lower = 0.8, upper = 0.95> pt1;
// real<lower = 365.25/4, upper = 365.25/2> tau;
// real<lower = 0.7, upper = 0.9> pst;
// real<lower = 0.4, upper = 0.6> RDT;
// real<lower = 0.1, upper = 0.15> ps;
// real<lower = 365.25/6, upper = 365.25/4> ss;
// real<lower = 365.25/7, upper = 365.25/3> omega;
// real<lower = 0.95, upper = 1.0> pt2;
// real<lower = 0.0005, upper = 0.002> trf1;
// real<lower = 0.0005, upper = 0.002> trf2;
// real<lower = 4000, upper = 6000> Kc;
// real<lower = 0.3, upper = 0.7> betaMH;
// real<lower = 0.3, upper = 0.7> betaHM;
// real<lower = 0.5, upper = 0.7> itn_eff;
// real<lower = 0.05, upper = 0.15> itn_aversion;
// real<lower = 0.3, upper = 0.5> itn_use;
// real<lower = 0.2, upper = 0.4> irs_eff;
// real<lower = 0.2, upper = 0.4> irs_aversion;
//standard deviations
real <lower=0> sigma_ICH; //Uncomplicated SD
real <lower=0> sigma_Sev; //Severe SD
}
transformed parameters {
real rains[T];
real temps[T];
real bites_rate[T];
real GammaH[T];
real chiM[T];
real muAM[T];
real PESA[T];
real muM[T];
real BT[T];
real upsilonEA[T];
real xiM[T];
real sigmaM[T];
real lambdaHM[T];
real lambdaMH[T];
real psi;
real mu;
for (t in 1:T) {
rains[t] = rainfall[t];
temps[t] = temperature[t];
bites_rate[t] = bites[t] + 9;
GammaH[t] = (1 - itn_cov[t] * itn_eff * (1 - itn_aversion) * itn_use) *
(1 - irs_cov[t] * irs_eff * (1 - irs_aversion));
chiM[t] = 30.44 * (-0.153 * pow(temps[t], 2) + 8.61 * temps[t] - 97.7);
muAM[t] = 30.44 * (1 / (8.560 + 20.654 * pow(1 + pow(temps[t] / 19.795, 6.827), -1)));
PESA[t] = (-0.00924 * pow(temps[t], 2) + 0.453 * temps[t] - 4.77);
muM[t] = 30.44 * (-log(-0.000828 * pow(temps[t], 2) + 0.0367 * temps[t] + 0.522));
BT[t] = chiM[t] / muM[t];
upsilonEA[t] = 30.44 * (1 / (-0.00094 * pow(temps[t], 2) + 0.049 * temps[t] - 0.552));
xiM[t] = BT[t] * PESA[t] / upsilonEA[t];
sigmaM[t] = (-0.00083 * pow(temps[t], 2) + 0.044 * temps[t] - 0.487);
lambdaHM[t] = ((bites_rate[t] * betaHM * (ICH + 0.5 * IAH + Sev +
0.5 * TACTH + 0.5 * TASH +
0.5 * TFACTH + 0.5 * TFASH)) / (SH + EH + ICH + IAH + Sev + TACTH + TASH + TFACTH + TFASH)) * GammaH[t];
lambdaMH[t] = ((bites_rate[t] * betaMH * IM) / (SH + EH + ICH + IAH + Sev + TACTH + TASH + TFACTH + TFASH)) * GammaH[t];
}
psi = pow((1 / lambdaMH[T] + 1 / gamma), -1);
mu = 1000 / 26.8;
}
model {
real dSH;
real dEH;
real dICH;
real dIAH;
real dSev;
real dTACTH;
real dTASH;
real dTFACTH;
real dTFASH;
real dAM;
real dSM;
real dEM;
real dIM;
for (t in 1:T) {
dSH = (SH + EH + ICH + IAH + Sev + TACTH + TASH + TFACTH + TFASH) / mu - m2 * lambdaMH[t] * SH + (1 - (ah1 + rs1)) * rho1 * TACTH + (1 - (ah2 + rs2)) * rho2 * TASH - SH / mu;
dEH = m2 * lambdaMH[t] * SH - pc2 * gamma * EH - (1 - pc2) * gamma * EH - EH / mu;
dICH = pc2 * gamma * EH + psi * IAH + q * IAH - pt1 * tau * pst * RDT * ICH - (1 - pst) * ps * ss * ICH - (1 - pst) * (1 - ps) * omega * ICH - ICH / mu;
dIAH = (1 - pc2) * gamma * EH + (1 - pst) * (1 - ps) * omega * ICH - psi * IAH - q * IAH - IAH / mu;
dSev = (1 - pst) * ps * ss * ICH - pt2 * ss * Sev - Sev / mu;
dTACTH = pt1 * tau * pst * RDT * ICH + trf1 * TACTH - (ah1 + rs1) * TACTH - (1 - (ah1 + rs1)) * rho1 * TACTH - TACTH / mu;
dTASH = pt2 * ss * Sev + trf2 * TASH - (ah2 + rs2) * TASH - (1 - (ah2 + rs2)) * rho2 * TASH + TFASH - TACTH / mu;
dTFACTH = (ah1 + rs1) * TACTH - trf1 * TFACTH - TFACTH / mu;
dTFASH = (ah2 + rs2) * TASH - TFASH - TFASH / mu;
dAM = chiM[t] * (1 - AM / (Kc * (rains[t] + 1))) * (SH + EH + ICH + IAH + Sev + TACTH + TASH + TFACTH + TFASH) - (xiM[t] + muAM[t]) * AM;
dSM = xiM[t] * AM - (lambdaHM[t] + muM[t]) * SM;
dEM = lambdaHM[t] * SM - (sigmaM[t] + muM[t]) * EM;
dIM = sigmaM[t] * EM - muM[t] * IM;
// Priors for new parameters
sigma_ICH ~ cauchy(0, 5);
sigma_Sev ~ cauchy(0, 5);
// Likelihood
//for (t in 1:T) {
observed_ICH[t] ~ normal(ICH, sigma_ICH);
observed_Sev[t] ~ normal(Sev, sigma_Sev);
}
}
generated quantities {
real ICH_sim[T];
real Sev_sim[T];
real log_lik[T];
for (t in 1:T) {
ICH_sim[t] = normal_rng(ICH, sigma_ICH);
Sev_sim[t] = normal_rng(Sev, sigma_Sev);
log_lik[t] = normal_lpdf(observed_ICH[t] | ICH, sigma_ICH) +
normal_lpdf(observed_Sev[t] | Sev, sigma_Sev);
}
}
stan Data
stan_data ← list(
T = length(UE$times),
times = UE$times,
rainfall=rainfall_fun(UE$times),
temperature=temperature_fun(UE$times),
bites=biting_fun(UE$times) + 9,
itn_cov=get_itn_coverage(UE$times),
irs_cov=get_irs_coverage(UE$times),
observed_ICH = UE$unc,
observed_Sev= UE$sev
)
Stan Fit
fit ← stan(file = “Model.stan”,
data = stan_data,
init = “random”,
iter = 5000,
warmup = 2500, # Half of total iterations
chains = 4,
cores = 4,
control = list(adapt_delta = 0.99, max_treedepth = 15)) # Adjust sampler parameters
posterior ← extract(fit)
simulated_ICH ← apply(posterior$ICH_sim, 2, mean)
simulated_Sev ← apply(posterior$Sev_sim, 2, mean)
Create a data frame for plotting
plot_data ← data.frame(
time = UE$times,
observed_ICH = observed_data$uncompmalaria,
simulated_ICH = simulated_ICH,
observed_Sev = observed_data$severemalaria,
simulated_Sev = simulated_Sev
)
Plot for ICH
ggplot(plot_data, aes(x = time)) +
geom_line(aes(y = observed_ICH/10000), linewidth = 0.5, color = “black”) +
geom_point(aes(y = observed_ICH/10000)) +
geom_line(aes(y = simulated_ICH/10000), color = “red”) +
labs(x = “Time”, y = “Cases ICH”) +
ggtitle("Combined Plot)Code / Program output
Plot for Sev
ggplot(plot_data, aes(x = time)) +
geom_line(aes(y = observed_Sev/10000), linewidth = 0.5, color = “black”) +
geom_line(aes(y = simulated_Sev/10000), color = “red”) +
labs(x = “Time”, y = “Cases Sev”) +
ggtitle(“Combined Plot”)