Programming a simple binomial GLM with logit link

Modeling
1

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
2

This should be C ~ binomial(N, inv_logit(logit_p));.

3

Thanks! I’ll try this and let you know if I have any further problems. As is so often the case, the answer is simpler than I thought.