Hello everyone, I am Ohene. I am new to Stan and R in general. I have only used both languages a few time but for simpler tasks. I am trying to estimate a linear stochastic response plateau with linear time trends to evaluate the effect of nitrogen application on wheat. The model follows the normal plateau model but has each coefficient (intercept, slope, and plateau) multiplied by the year variable in addition to been multiplied by the nitrogen variable where applicable and year random effects. Unfortunately, I cannot share the data due to confidential reasons, but the three major variables will be the “year” (1970-2010), “yield”, and nitrogen (0,50,200,150,100). I run the model in R with the code below. I get the error
“Error in stanc(file = file, model_code = model_code, model_name = model_name, :
parser failed badly”.
I am trying to understand the error but also want to know if I did anything wrong in the model build up in stan or any error i could have made. I would be grateful for your explanations and help. (PS, I apologize if i pasted in the wrong forum)
random_effects <- '
data {
int N; // number of obs
int M; // number of groups (years)
int g[N]; // map obs to groups
vector[N] y; //the response variable
vector[N]x; // predictor or nitrogen
vector [N]Q; // year variable
}
parameters {
real a[M];
real<lower=0> b0; // intercept
real<lower=0> b0_time; //intercept linear trend
real<lower=0> b1 ; //slope
real<lower=0> b1_time; //slope linear trend
real<lower=0> plateau; // plateau level
real<lower=0> plateau_time; // plateau level time trend
real<lower=0> plateau_error; // plateau
real<lower=0> year_error;
real<lower=0> s2u;
real<lower=0> s2v;
real<lower=0> sigmae;
}
model {
# Priors
b0 ~ normal(147, 1e6);
b0_time ~normal(0, 1e6);
b1 ~ normal(0.45, 1e6);
b1_time ~ normal(0, 1e6);
plateau ~ normal(183.55, 1e6);
plateau_time ~ normal(0, 1e6);
plateau_error ~ normal(0,s2u);
year_error ~ normal(0,s2v);
s2u ~ gamma(0.001,10000);
s2v ~ gamma(0.001,10000);
sigmae ~ gamma(0.001,10000);
a ~ normal(0,1);
for (j in 1:N) {
response[j]=fmin((b0 + b0_time*Q[j] + (b1 + b1_time*Q[j])*x[j]), plateau + plateau_time*Q[j] + plateau_error[g[j]]) + year_error*a[g[j]];
}
for(n in 1:N) {
y[n] ~ normal( response[n],sigmae); //R uses standard deviation
}
}
'
library(parallel)
getOption("mc.cores", 1L) #Does nothing here since cores is specified below
library(rstan)
library(readxl)
library(dplyr)
setwd("C:\\Users\\ohene\\OneDrive-KNUST \\Documents\\Wheat Response to Nitrogen\\Data and Codes\\Data")
data<-read_excel("Data.xlsx")
data <- as.data.frame(data)
nobs <- nrow(data)
years <- length(unique(data[,"YR"]))
xx <- data$N
gg <- group_indices(data,YR) #creates an index for year
y <- data$BUAC
qq <- data$YR
xx <- cbind(xx) #intercepts are added in the code
qq <-cbind(qq)
yield_data <- list(N=NROW(y),M=years,K=1,y=y,x=xx,Q=qq,g=gg)
hfit <- stan(model_code=random_effects, model_name="lrsp",data=yield_data, iter=3000, warmup=1000, chains=8, cores=8, seed = 29473, control = list(adapt_delta = 0.80, max_treedepth = 15)).