I’m working my way through some of the codes in Kery&Schub’s Bayesian Pop analysis Using winbugs as posted on Stan by Trevor Branch. The example is from Chapter 3, section 3.5.2 analyzing real data (falcons attached). I ran the code for GLM_Binomial.stan Trevor provided (with a little editing) just fine. This states the model likelihood as C ~ binomial_logit(N, logit_p). Although I understand the efficiency of this coding, I wanted to also code it directly using binomial(N,p). But however, I have tried it I get error messages upon trying to compile in Rstan, that seem to have to do with declaring p first and logt_p as a derived parameter, rather than the other way around (using binomial_logit). I’d just like to know if it is possible to code it the other way around declaring p first and using binomial(N, p) while, of course, defining the logit link. The data, a test .stan file, and the associated .R file are attached.

Thanks for any help, Nick

data {

int<lower=0> nyears; // Number of Years

int<lower=0> C[nyears]; // Counts

int<lower=0> N[nyears]; // Binomial Totals

vector[nyears] yyear; // Year covariates

}

transformed data {

vector[nyears] year_squared;

year_squared = yyear .* yyear;

}

parameters {

real alpha;

real beta1;

real beta2;

real<lower=0, upper=1> p[nyears];

}

transformed parameters {

real logit_p[nyears];

// Linear predictor

logit_p = logit§;

logit_p = alpha

+ beta1 * yyear

+ beta2 * year_squared;

}

model {

// Priors

alpha ~ normal(0, 100);

beta1 ~ normal(0, 100);

beta2 ~ normal(0, 100);

// Likelohood

C ~ binomial(N, p);

// C ~ binomial_logit(N, logit_p); original likelihood

}

//generated quantities {

// real<lower=0,upper=1> p[nyears];

// for (i in 1:nyears)

// p[i] = inv_logit(logit_p[i]);

//}

## 3. Introduction to the generalized linear model (GLM): The simplest

## model for count data

## 3.5. Binomial GLM for modeling bounded counts or proportions

## 3.5.1. Generation and analysis of simulated data

library(rstan)

rstan_options(auto_write = TRUE)

## Read data

Birds <- read.table(“falcons.txt”, header = TRUE)

C<-Birds$RPairs

N<-Birds$Count

Year<-Birds$Year

nyears <-length©

yyear<-(Year-1985)/20

stan_data<- list(C, N, nyears, yyear)

## Call Stan from R

fit <- stan(“GLM_Binom_Test.stan”, data = stan_data,init = list(list(alpha = 1.0, beta1 = 1.0, beta2 = 1.0)), iter = 100, chains = 1)

## Summarize posteriors

print(fit)

Year | Count | RPairs | Eyasses |
---|---|---|---|

1964 | 34 | 19 | 27 |

1965 | 45 | 18 | 30 |

1966 | 39 | 22 | 36 |

1967 | 36 | 19 | 37 |

1968 | 20 | 8 | 19 |

1969 | 18 | 7 | 10 |

1970 | 21 | 9 | 9 |

1971 | 20 | 11 | 17 |

1972 | 17 | 10 | 15 |

1973 | 20 | 12 | 23 |

1974 | 23 | 15 | 35 |

1975 | 27 | 19 | 40 |

1976 | 29 | 21 | 42 |

1977 | 39 | 27 | 53 |

1978 | 43 | 32 | 65 |

1979 | 47 | 35 | 77 |

1980 | 57 | 42 | 97 |

1981 | 56 | 36 | 81 |

1982 | 63 | 52 | 112 |

1983 | 70 | 44 | 85 |

1984 | 82 | 58 | 121 |

1985 | 90 | 54 | 117 |

1986 | 95 | 59 | 122 |

1987 | 107 | 72 | 154 |

1988 | 114 | 80 | 141 |

1989 | 121 | 81 | 173 |

1990 | 132 | 95 | 205 |

1991 | 134 | 101 | 209 |

1992 | 143 | 92 | 176 |

1993 | 150 | 114 | 238 |

1994 | 154 | 74 | 143 |

1995 | 172 | 79 | 121 |

1996 | 163 | 115 | 233 |

1997 | 164 | 102 | 225 |

1998 | 161 | 102 | 200 |

1999 | 137 | 76 | 112 |

2000 | 157 | 104 | 240 |

2001 | 153 | 64 | 130 |

2002 | 191 | 107 | 178 |

2003 | 193 | 104 | 197 |