Workshop 10.6b: Analysis of count data (Bayesian) Murray Logan - PDF document

-1- Workshop 10.6b: Analysis of count data (Bayesian) Murray Logan September 13, 2016 Table of contents 1 Poisson regression 1 2 Helper Functions 6 3 Worked Examples 7 1. Poisson regression 1.1. Poisson regression Probability density function Cumulative density function λ = 25 λ = 15 λ = 3 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 p ( Y i ) = e − λ λ x x ! log ( µ ) = β 0 + β 1 x i + ... + β p x p 1.2. Dispersion Spread assumed to be equal to mean. ( φ = 1 )

-2- Probability density function Cumulative density function λ = 25 λ = 15 λ = 3 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 1.3. Dispersion 1.3.1. Over-dispersion Sample more varied than expected from its mean • variability due to other unmeasured ( latent ) influences – quasi-poisson model – negative binomial – observation level random effect 1.4. Dispersion 1.4.1. Over-dispersion Sample more varied than expected from its mean • variability due to other unmeasured ( latent ) influences • clumpiness – negative binomial model • due to more zeros than expected – zero-inflated model 1.5. Residuals • difficult to interpret

-3- 1.6. Simulated Residuals • residuals uniform when model good Raw Pearson ● 4 3 ● ● ● ● ● ● ● ● ● ● ● ● ● 3 2 ● ● ● ● ● 2 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 0 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1 −1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2 height! height! 1.68 1.70 1.72 1.74 1.76 1.78 1.68 1.70 1.72 1.74 1.76 1.7. Simulated Residuals • residuals uniform when model good

-4- Raw Pearson ● 4 3 ● ● ● ● ● ● ● ● ● ● ● ● ● 3 2 ● ● ● ● ● 2 ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 0 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1 −1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2 height! height! 1.68 1.70 1.72 1.74 1.76 1.78 1.68 1.70 1.72 1.74 1.76 Raw Pearson ● 4 3 ● ● ● ● ● ● ● ● ● ● ● ● ● 3 2 ● ● ● ● ● 2 ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 0 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1 −1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2 height! height! 1.68 1.70 1.72 1.74 1.76 1.78 1.68 1.70 1.72 1.74 1.76 1.8. Simulated Residuals • residuals uniform when model good

-5- Raw Pearson ● 4 3 ● ● ● ● ● ● ● ● ● ● ● ● ● 3 2 ● ● ● ● ● 2 ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 0 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1 −1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2 height! height! 1.68 1.70 1.72 1.74 1.76 1.78 1.68 1.70 1.72 1.74 1.76 Raw Pearson ● 4 3 ● ● ● ● ● ● ● ● ● ● ● ● ● 3 2 ● ● ● ● ● 2 ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 0 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1 −1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2 height! height! 1.68 1.70 1.72 1.74 1.76 1.78 1.68 1.70 1.72 1.74 1.76 1.9. Simulated Residuals • residuals uniform when model good

-6- Raw Pearson ● 4 3 ● ● ● ● ● ● ● ● ● ● ● ● ● 3 2 ● ● ● ● ● 2 ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 0 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1 −1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2 height! height! 1.68 1.70 1.72 1.74 1.76 1.78 1.68 1.70 1.72 1.74 1.76 Raw Pearson ● 4 3 ● ● ● ● ● ● ● ● ● ● ● ● ● 3 2 ● ● ● ● ● 2 ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 0 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1 −1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −2 height! height! 1.68 1.70 1.72 1.74 1.76 1.78 1.68 1.70 1.72 1.74 1.76 2. Helper Functions

-7- 2.1. Simulated Residuals > simRes <- function(lambda, Y,n=250, plot=T, family='negbin', size=NULL,theta=NULL) { + require(gap) + N = length(Y) + sim = switch(family, + 'poisson' = matrix(rpois(n*N,apply(lambda,2,mean)),ncol=N, byrow=TRUE), + 'negbin' = matrix(MASS:::rnegbin(n*N,apply(lambda,2,mean),size),ncol=N, byrow=TRUE), + 'zip' = matrix(gamlss.dist:::rZIP(n*N,apply(lambda,2,mean),theta),ncol=N, byrow=TRUE), + 'zinb' = matrix(gamlss.dist:::rZINBI(n*N,apply(lambda,2,mean),sigma=theta,nu=size),ncol=N, + byrow=TRUE) + ) + a = apply(sim + runif(n,-0.5,0.5),2,ecdf) + resid<-NULL + for (i in 1:length(Y)) resid<-c(resid,a[[i]](Y[i] + runif(1 ,-0.5,0.5))) + if (plot==T) { + par(mfrow=c(1,2)) + gap::qqunif(resid,pch = 2, bty = "n", + logscale = F, col = "black", cex = 0.6, main = "QQ plot residuals", + cex.main = 1, las=1) + plot(resid~apply(lambda,2,mean), xlab='Predicted value', ylab='Standardized residual', las=1) + } + resid + } 2.2. ggplot goodness > library(ggplot2) > theme_goodness <- theme_classic() + + theme(axis.line.x=element_line(), axis.line.y=element_line()) 2.3. MCMC summaries > MCMCsum <- function(x) { + data.frame(Mean=mean(x), median=median(x), HPDinterval(as.mcmc(x))) + } > MCMCsummaries <- function(fit) { + plyr:::adply(fit, 2, MCMCsum) + } 3. Worked Examples 3.1. Poisson t-test Format of limpets.csv data files Count Shore 1 sheltered 3 sheltered 2 sheltered 1 sheltered 4 sheltered . . .