Gaussian vs Poisson Regression

MC simulation and the boostrap

1 The data

This data is taken from Section 8.8 (Modeling Count Data) from Jim Albert and Jingchen Hu’s 2020 book Probability and Bayesian Modeling (CRC Press, Taylor & Francis Group). Suppose one is interested in monitoring the popularity of a particular blog about baseball analytics. The data in y (below) displays the number of visitors viewing this blog for 28 days during June of 2018.

y = c(95,81,85,100,111,130,113,92,65,78,96,118,120,104,91,91,79,106,91,114,110,98,61,84,96,126,119,90)

n = length(y)

days = c("fri","sat","sun","mon","tue","wed","thu")

t = rep(1:7,4)

# Exploratory data analysis

par(mfrow=c(2,2))
hist(y,xlab="",main="Number of visitors")
boxplot(y,ylab="Number of visitors")

ts.plot(y,type="b",xlab="Days",ylab="Number of visitors")
abline(v=7.5,lty=2)
abline(v=14.5,lty=2)
abline(v=21.5,lty=2)
 
plot(t,y,xlab="Day of the week",ylab="Number of visitors",axes=FALSE)
axis(2);axis(1,at=1:7,lab=days)

2 Linear regression vs Poisson regression

weekend=rep(0,n)
weekend[t<=3]=1

c(mean(y[weekend==0]),sqrt(var((y[weekend==0]))))

## [1] 109.00000  12.15456

c(mean(y[weekend==1]),sqrt(var((y[weekend==1]))))

## [1] 83.33333 11.42034

par(mfrow=c(1,2))
boxplot(y,ylab="Number of visitors")
boxplot(y~weekend,xlab="Type of day",names=c("Weekday","Weekend"),ylab="Number of visitors")

2.1 Gaussian regression

summary(lm(y~weekend))

## 
## Call:
## lm(formula = y ~ weekend)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.333  -6.250   1.333   8.750  21.000 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  109.000      2.962  36.795  < 2e-16 ***
## weekend      -25.667      4.525  -5.672 5.76e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.85 on 26 degrees of freedom
## Multiple R-squared:  0.5531, Adjusted R-squared:  0.5359 
## F-statistic: 32.17 on 1 and 26 DF,  p-value: 5.76e-06

2.2 Poisson regression

summary(glm(y~weekend,family="poisson"))

## 
## Call:
## glm(formula = y ~ weekend, family = "poisson")
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.69135    0.02395 195.916  < 2e-16 ***
## weekend     -0.26850    0.03967  -6.769  1.3e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 85.141  on 27  degrees of freedom
## Residual deviance: 38.324  on 26  degrees of freedom
## AIC: 221.76
## 
## Number of Fisher Scoring iterations: 4

tab = rbind(c(109.000,109.000-25.667),
c(exp(4.69135),exp(4.69135-0.26850)))

rownames(tab) = c("Gaussian regression","Poisson regression")
colnames(tab) = c("Intercept","Slope")
tab

##                     Intercept    Slope
## Gaussian regression  109.0000 83.33300
## Poisson regression   109.0002 83.33345

ma = 4.7;sda=0.1
mb = -0.3;sdb=0.1
N = 100
alpha = seq(4.5,4.9,length=N)
beta = seq(-0.5,0,length=N)
loglike = matrix(0,N,N)
prior = matrix(0,N,N)
posterior = matrix(0,N,N)
for (i in 1:N)
  for (j in 1:N){
    lambda = exp(alpha[i]+beta[j]*weekend)
    loglike[i,j] = sum(dpois(y,lambda,log=TRUE))
    prior[i,j] = dnorm(alpha[i],ma,sda)*dnorm(beta[j],mb,sdb)
  }
like = exp(loglike)
post = prior*like

par(mfrow=c(1,1))
image(alpha,beta,prior,xlab=expression(alpha),ylab=expression(beta))
contour(alpha,beta,like,col=4,drawlabels=FALSE,add=TRUE,xlab=expression(alpha),ylab=expression(beta),lwd=2)
contour(alpha,beta,post,col=3,drawlabels=FALSE,add=TRUE,xlab=expression(alpha),ylab=expression(beta),lwd=2)

legend("topright",legend=c("Prior","Likelihood","Posterior"),col=c(2,4,3),lty=1,lwd=3)

3 SIR-based posterior inference

M = 10000
alphas = rnorm(M,ma,sda)
betas  = rnorm(M,mb,sdb)
plot(alphas,betas,xlab=expression(alpha),ylab=expression(beta))
contour(alpha,beta,prior,add=TRUE,drawlabels=FALSE,lwd=2,col=2)

w = rep(0,M)
for (i in 1:M)
  w[i] = sum(dpois(y,exp(alphas[i]+betas[i]*weekend),log=TRUE))
w = exp(w-max(w))
ind = sample(1:M,size=M,replace=TRUE,prob=w)

alphas1 = alphas[ind]
betas1 = betas[ind]
plot(alphas1,betas1)

image(alpha,beta,post,xlab=expression(alpha),ylab=expression(beta))
contour(alpha,beta,post,add=TRUE,drawlabels=FALSE,lwd=2,col=2)

par(mfrow=c(1,1))
plot(density(exp(alphas1)),xlim=c(70,120),ylim=c(0,0.2),xlab="Rate of visitors",ylab="Density",main="",lwd=2)
lines(density(exp(alphas1+betas1)),col=2,lwd=2)
legend("topright",legend=c("Weekdays","Weekends"),col=1:2,lty=1,lwd=2)