// generated with brms 2.10.0 functions { /* Temperature multiplier */ real temFun(real tem) { real Tmin; real Topt; real Tmax; real alpha; Tmin = 0; Topt = 30; Tmax = 45; alpha = log(2)/log((Tmax - Tmin)/(Topt - Tmin)); if (tem <= Tmin || tem >= Tmax) return 0; else return (2 * ((tem - Tmin) ^ alpha) * ((Topt - Tmin) ^ alpha) - ((tem - Tmin) ^ (2 * alpha)))/ ((Topt - Tmin) ^ (2 * alpha)); } /* Photoperiod multiplier */ real phoFun(real pho, real Pcrit, real Psen) { real x; if(pho <= Pcrit) x = 1; else x = 1 - Psen * (pho - Pcrit); if (x < 0) x = 0; return x; } } data { int N; // number of observations vector[N] Y; // response variable int Nmi; // number of missings int Jmi[Nmi]; // positions of missings int K_dOptFiFl; // number of population-level effects matrix[N, K_dOptFiFl] X_dOptFiFl; // population-level design matrix int K_dOptFlSd; // number of population-level effects matrix[N, K_dOptFlSd] X_dOptFlSd; // population-level design matrix int K_dOptSdPm; // number of population-level effects matrix[N, K_dOptSdPm] X_dOptSdPm; // population-level design matrix int K_Pcrit; // number of population-level effects matrix[N, K_Pcrit] X_Pcrit; // population-level design matrix int K_Psen; // number of population-level effects matrix[N, K_Psen] X_Psen; // population-level design matrix int KC; // number of covariates matrix[N, KC] C; // matrix of covariates int prior_only; // should the likelihood be ignored? int time[N]; // time index } transformed data { // extract covariate vectors for faster indexing vector[N] C_1 = C[, 1]; vector[N] C_2 = C[, 2]; real dOptEmJuv; real dOptJuvFi; dOptEmJuv = 3.7; dOptJuvFi = 3.7; } parameters { vector[Nmi] Ymi; // estimated missings vector[K_dOptFiFl] b_dOptFiFl; // population-level effects vector[K_dOptFlSd] b_dOptFlSd; // population-level effects vector[K_dOptSdPm] b_dOptSdPm; // population-level effects vector[K_Pcrit] b_Pcrit; // population-level effects vector[K_Psen] b_Psen; // population-level effects real sigma; // residual SD real init; // initial value of mu } transformed parameters { // vector combining observed and missing responses vector[N] Yl = Y; // initialize linear predictor term vector[N] nlp_dOptFiFl = X_dOptFiFl * b_dOptFiFl; // initialize linear predictor term vector[N] nlp_dOptFlSd = X_dOptFlSd * b_dOptFlSd; // initialize linear predictor term vector[N] nlp_dOptSdPm = X_dOptSdPm * b_dOptSdPm; // initialize linear predictor term vector[N] nlp_Pcrit = X_Pcrit * b_Pcrit; // initialize linear predictor term vector[N] nlp_Psen = X_Psen * b_Psen; // initialize non-linear predictor term vector[N] mu; Yl[Jmi] = Ymi; mu[1] = init; for (n in 2:N) { // compute non-linear predictor values if (time[n] == 1) { mu[n] = init; } else if (time[n] > 1 && mu[n] <= 1) { mu[n] = mu[n - 1] + 1 / dOptEmJuv * temFun(C_1[n]); } else if (time[n] > 1 && mu[n] > 1 && mu[n] <= 2) { mu[n] = mu[n - 1] + 1 / dOptJuvFi * temFun(C_1[n]) * phoFun(C_2[n] , nlp_Pcrit[n], nlp_Psen[n]); } else if (time[n] > 1 && mu[n] > 2 && mu[n] <= 3) { mu[n] = mu[n - 1] + 1 / nlp_dOptFiFl[n] * temFun(C_1[n]) * phoFun(C_2[n] , nlp_Pcrit[n], nlp_Psen[n]); } else if (time[n] > 1 && mu[n] > 3 && mu[n] <= 4) { mu[n] = mu[n - 1] + 1 / nlp_dOptFlSd[n] * temFun(C_1[n]) * phoFun(C_2[n] , nlp_Pcrit[n], nlp_Psen[n]); } else { mu[n] = mu[n - 1] + 1 / nlp_dOptSdPm[n] * temFun(C_1[n]) * phoFun(C_2[n] , nlp_Pcrit[n], nlp_Psen[n]); } } } model { // priors including all constants target += normal_lpdf(init | 0, 0.01); target += normal_lpdf(b_dOptFiFl | 15, 6); target += normal_lpdf(b_dOptFlSd | 16, 7); target += normal_lpdf(b_dOptSdPm | 40, 9); target += normal_lpdf(b_Pcrit | 13, 1); target += normal_lpdf(b_Psen | 0.29, 0.05); target += normal_lpdf(sigma | 0.0, 0.1); // likelihood including all constants target += normal_lpdf(Yl | mu, sigma); } generated quantities { }