# GLM and implementation of Poisson regression model in R by hand

first of, this is not my school exercise but a given example that I'd like to convert from Stan to my own code. I am very much a pragmatic learner so doing this helps me a lot to visualize the problem.

The problem is the following:

#The data give numbers of fatal accidents on scheduled airline flights per year over a ten-year
#period. Assume that the number of fatal accidents in year t follows a Poisson distribution with
#mean theta where log(theta)=a+bt.

#Year      1976 1977 1978 1979 1980 1981 1982 1983 1984 1985
#accidents  24  25    31   31   22   21   26   20   16   22

#theta[i]  # true' mean number of fatal accidents
y = c(24, 25, 31, 31, 22, 21, 26, 20, 16, 22)                      # number of fatal accidents
t = c(1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985) # year

#Using noninformative priors on a and b, obtain the posterior distributions of a and b. Plot
#the approximate probability density for the expected number of fatal accidents in year 1988.
#Obtain the predictive interval for the the number of fatal accidents in 1988.


Which is run with the following Stan model:

//Poisson_regression_model
data{
int N; //the number of observations
real x[N]; //years, real because of mean
int y[N]; // accidents
}
parameters{
real a_star;
real b;
}
transformed parameters{
vector[N] log_theta;
vector[N] theta;
real a;
a=a_star-b*mean(x);

for( i in 1 : N ) {
log_theta[i] = a_star+b*(x[i]-mean(x));
}

theta=exp(log_theta);
}
model{
//priors
a_star ~ normal(0, 10^3);
b ~ normal(0, 10^3);

//likelihood
for( i in 1 : N ) {
//see Stan reference manual page 520
y[i] ~ poisson_log(log_theta[i]);
}
}


Full code and which I'd like to convert into my own, very naive implementation. So this is how far I got but I must have some mistake calculating my likelihood/log_theta.

a_param_space <- seq(0.1, 20, length=100)
b_param_space <- seq(0.1, 20, length=100)
a_prior <- function(x) dnorm(x, 0, 1e3)
b_prior <- function(x) dnorm(x, 0, 1e3)

theta_likelihood <- function(x, theta) dpois(x, theta, log = T)

y = c(24, 25, 31, 31, 22, 21, 26, 20, 16, 22)
t = c(1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985)

posterior <- function() {
plen <- length(a_param_space)
res <- matrix(rep(0, plen*2), ncol=2)
for(i in 1:plen) {
for(j in 1:plen) {
a <- a_param_space[i]
b <- b_param_space[j]
log_theta <- a + b*(t - mean(t))
print(log_theta)
ab_likelihood <- theta_likelihood(y, log_theta)
print(ab_likelihood)
ab_prior_likelihood <- sum(a_prior(a) + b_prior(b) + ab_likelihood)
print(ab_prior_likelihood)
res[i, j] = ab_prior_likelihood
}
}
a_sum <- sum(res[,1])
b_sum <- sum(res[,2])
normalized <- matrix(rep(0, plen*2), ncol=2)
normalized[,1] <- res[,1]/a_sum
normalized[,2] <- res[,2]/b_sum
return(normalized)
}


Also I'm not sure have I named my variables correctly, it's very hard to me to convert the mathematical formula into programming terms. Say as another example I have a likelihood function $$y_i|\theta_i \sim Bin(n_i, \theta_i)$$ where $$logit(\theta_i)=\alpha + \beta x_i$$ with the link function $$logit(\theta_i)=log(\theta_i/(1-\theta_i)$$.

So this would mean I have to pass as a parameter to binomial density the values y-vector, n-vector and theta, where theta is computed from two parameters (with some values from parameter space), alpha and beta. Those two are then computed from the linear regression model: alpha + beta * x-vector which are inputted to the logit function (alpha + beta * x-vector)/(1-alpha + beta * x-vector) to get as a result a vector of log-odds? Am I almost correct?

Thank you in advance for answering, I am very much a caricature programmer who knows how to code but can't understand math. No need to rub it in.

• Hi TeemuK - have you looked into the use of the rstanarm package in R? The reason I ask is that you might find it simpler to specify your model in R using the built functions related to setting priors for your model parameters and for specifying your likelihood function in that syntax format first. You can have the function write the stan code to file for inspection. Also you are using pretty diffuse priors - the Stan manual recommends weakly informative priors whenever possible. – Matt Barstead Dec 8 '18 at 16:39
• You might also find it easier to specify the linear predictor first in your transformed parameters statement, then exponentiate, and then model that term as Poisson-distributed. See this link for an example: rpubs.com/kaz_yos/stan-pois1 – Matt Barstead Dec 8 '18 at 17:11

Well I did manage to create the Poisson regression using MAP estimation (not MLE :)). Although I'm sure there is still some parts of it wrong, I'm happy that I managed to pull it off. It has pretty funky looking posterior distribution. However, I didn't figure out how to do the binomial regression. Dunno will I bother looking into it.

a_param_space <- seq(2, 4, length=100)
b_param_space <- seq(-1, 1, length=100)
a_prior <- function(x) dnorm(x, 0, 1e3)
b_prior <- function(x) dnorm(x, 0, 1e3)

poisson_log_ll <- function(y, x, a, b) {
theta = exp(a + b*x) # link function
return(sum(dpois(y, theta), log=TRUE))
}

y = c(24, 25, 31, 31, 22, 21, 26, 20, 16, 22)
t = c(1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985)
tt = t - 1976 # Normalize to 0

posterior <- function() {
plen <- length(a_param_space)
res <- matrix(rep(0, plen*plen), ncol=plen)
for(i in 1:plen) {
for(j in 1:plen) {
a <- a_param_space[i]
b <- b_param_space[j]
ab_likelihood <- poisson_log_ll(y, tt, a, b)
ab_prior_likelihood <- sum(a_prior(a) + b_prior(b) + ab_likelihood)
res[i, j] = ab_prior_likelihood
}
}
normalized <- matrix(rep(0, plen*plen), ncol=plen)
normalized <- res/sum(res)
return(normalized)
}

map <- function(post) {
res = which(post==max(post), arr.ind = TRUE)
return(c(a_param_space[res], b_param_space[res]))
}

pois_post <- posterior()
print(map(pois_post))
print(glm(y ~ tt, family=poisson)\$coefficients)

contour(y=a_param_space, x=b_param_space, z=pois_post)
image(y=a_param_space, x=b_param_space, z=pois_post)
persp(y=a_param_space, x=b_param_space, z=pois_post, theta = 25, phi = 35)
persp(y=a_param_space, x=b_param_space, z=pois_post, theta = 0, phi = 90)
` 