“Precision (certainty) is important in making decisions” is a statement that I have little uncertainty that even the most disagreeable of people will accept. But how often do you check that your estimated uncertainty is consistent with your observed uncertainty?

A common barrier is that for GLMs it may not be clear what the model predicted variance should be. Which is why as always your selfless author has taken it upon him self to explain (i.e., I talk too much about things that no one else cares about for the gratification of my ego).

For example Fig. 1 shows histograms of two sets of quiz scores. The means of both sets are indistinguishable from their true mean of 5, however the variances are 21.30 and 2.72 respectively. The variance of a binomial is ; for and the variance should be 2.5, thus the first has a much higher observed variance then theoretical variance.

The ratio of the theoretical variance to the observed variance is called dispersion (). If then the model is over-dispersed; if then the model is under-dispersed. This post will consider only over-dispersion. Over-dispersion is evidence that there are missing explanatory variables.

In the quiz example there are clearly two groups in class 1 but they are modeled as only one group. Thus the best way to model over-dispersed data is to find the missing underling predictor however in general this might not be possible (or desirable if a simple uni-variate relationship is sought). The non-Bayesian way (i.e., inferior way) to handle dispersion is to use a quasi-likelihood GLM (Prof. Altman class is where I learned about this in my under grad and I still refer to her excellent notes; If you are reading this Dr Altman please convert your notes into a text book). To summarize this method, the log likelihood is replaced with a score function () that has the same properties (, , and ) used to derive the asymptotic results in exponential family GLM modeling. Here is the number of observations, is the model parameter, and is the data uncertainty. The parameters are then estimated optimizing instead of likelihood.

This method has the advantage that it can be implement in r with only a few extra key strokes:

“glm(cbind(Correct,Wrong )~factor(Class), data=data.example, family=binomial(link=”probit”) )”

is the normal way and

“glm(cbind(Correct,Wrong )~factor(Class), data=data.example, family=**quasi**binomial(link=”probit”) )”

is the quasi likelihood way. The down side is that quasi likelihood is not likelihood and consequently can not be used in many Bayesian model selection schemes (e.g., there is no evidence estimate or even BIC for quasi-likelihood models).

The Bayesian way (i.e., the correct way) is to use a prior distribution with non stated hyper parameters () to give the model the flexibility to account for dispersion. The “recipe” is to make . For the quiz example this means using a binomial distribution for and a beta distribution for the result of the integral is then the beta binomial distribution, though if you don’t want to take wiki’s word for it, this is a fun integral to work though in a pretentious coffee shop on a rain day. Particularly if you are trying to pull an attractive 3rd math student; the functions give it a pleasantly complex appearance. In other words, the beta binomial distribution is the product of a binomial distribution with a beta prior having the probability of a success integrated out. Conceptually this means that the probability of a success is random from a beta distribution; unknown predictors are then absorbed into this distribution.

The pdf of beta binomial is where is the probability of a success and is precision. For a GLM it would be normal to define (the function would typically be either a probit or logit) and estimate on its own. This is analogous to estimating the mean and sd for a linear regression.

The appended code gives an example of estimating the probability of success answering a question for both class with 4 additional confounding parameters. A classic binomial GLM is also run. For the GLM, most simulated data sets (about 80% but I have not conducted an exhaustive simulation) find at least one confounding parameters to be a significant predictor of probability of success. The GLM summary for one run is

“Call:

glm(formula = cbind(k, n – k) ~ ., family = binomial(link = “logit”),

data = X)

Deviance Residuals:

Min 1Q Median 3Q Max

-4.0594 -2.1550 -0.1402 2.1930 4.4550

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 0.137667 0.126606 1.087 0.276875

V1 -0.235301 0.181677 -1.295 0.195264

V2 -0.177382 0.085473 -2.075 0.037959 *

V3 0.078763 0.090477 0.871 0.384009

V4 0.351333 0.102759 3.419 0.000629 ***

V5 0.001959 0.092348 0.021 0.983073

—

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 381.54 on 59 degrees of freedom

Residual deviance: 366.40 on 54 degrees of freedom

AIC: 478.64

Number of Fisher Scoring iterations: 4″

The result of the Beta Binomial model is

The Beta binomial model correctly does not find any of the parameters to be significantly different from zero.

Code:

g.beta.binomial <- function(x, n, mu=0.5, phy=1000, Log.out =T )

{

#print(is.numeric(n))

#print(is.numeric(x))

#print(is.numeric(phy))

#print(is.numeric(mu))

#temp=data.frame(A1 = x+phy*mu, B1 = n-x+phy*(1-mu), A2 = phy*mu, B2= phy*(1-mu))

#print(temp)

#LP = lbeta(temp$A1, temp$B1)-lbeta( phy*mu, phy*(1-mu))+lchoose(n, x)

LP = lbeta(x+phy*mu, n-x+phy*(1-mu))-lbeta( phy*mu, phy*(1-mu))+lchoose(n, x)

if(Log.out == F) {LP = exp(LP)}

LP[is.nan(LP)] = -300 # get NaNs out

LP = sum(LP)

#print(LP)

return(LP)

}

g.LL <- function(M, K, n, X, npar, ndat) # the log likelihood

{

eta = 0*(1:ndat) + M[2] # the systimatic part

#print(is.numeric(eta))

for (i in 1:npar)

{

#print(is.numeric(M[2+i]))

#print(is.numeric(X[,i]))

eta = eta + M[2+i]*X[,i]

#print(is.numeric(eta))

}

Mu = 1/(1+exp(-eta)) # the expected value

Phy = M[npar+3]

#print(K)

#print(n)

LL = g.beta.binomial(K, n, mu=Mu, phy=Phy, Log.out =T )

return(LL)

}

# does not do rotation

g.BBstep <- function(M, k., n., X., Q.sd., npar., ndat.)

{

# int step

new = M

new[2] = new[2] + Q.sd.[1]*rnorm(1)

new[1] = g.LL(new, k., n., X., npar., ndat.)

r = runif(1)

A = exp(new[1]-M[1])

if (is.na(A)) {A= 0}

if (A>=r)

{

M = new

}

# the model pars

for (i in 1:npar.)

{

j = i + 2

new = M

new[j] = new[j] + Q.sd.[i+1]*rnorm(1)

new[1] = g.LL(new, k., n., X., npar., ndat.)

r = runif(1)

A = exp(new[1]-M[1])

if (is.na(A)) {A= 0}

if (A>=r)

{

M = new

}

}

# int step

new = M

new[npar.+3] = new[npar.+3] + Q.sd.[npar.+2]*rnorm(1)

if(new[npar.+3] > 0)

{new[1] = g.LL(new, k., n., X., npar., ndat.)

r = runif(1)

A = exp(new[1]-M[1])

if (is.na(A)) {A= 0}

if (A>=r)

{

M = new

}

}

return(M)

}

g.BB.glm <- function(k, n, X, NSTEPS = 50000, nsave = 1000)

{

npar = ncol(X) # get the number of pars

ndat = length(k)

parnames = colnames(X) # par names

old = (1:(npar+3))*0 # LL, intersept, pars, dispersion

# this gets the starting model and step size

MOD.start = glm(cbind(k, n-k)~., data=X, family = binomial(link = “logit”))

print(summary(MOD.start))

old[1] = -Inf

old[2:(npar+2)] = MOD.start$coefficients # staring model

old[npar+3] = 100*sum(n) # starting dispersion

Q.sd = sqrt(diag(vcov(MOD.start))) # step size

Q.sd = c(Q.sd, sum(n))

print(Q.sd)

# fixed length thinning

par.buf = matrix(0, nrow=nsave, ncol=(npar+1+1+1+1)) # the saved models pars

# this part samples and saves the models

print(date())

burn= 10000

for (i in 1:(nsave+burn)) # loop over sampling first

{

old = g.BBstep(old, k, n, X, Q.sd.=Q.sd, npar, ndat) # new model (1 swweep) with LL

#print(old)

if (i >= burn +1)

{

j = i – burn

par.buf.temp = c(j, old)

#print(par.buf.temp)

#print(par.buf[i,])

par.buf[j,] = par.buf.temp

#print(par.buf[i,])

}

}

print(date())

for (i in 3:(ncol(par.buf)))

{

Q.sd[i-2] = sd(par.buf[,i])

}

Q.sd[length(Q.sd)] = Q.sd[length(Q.sd)]

print(cor(par.buf[,3:ncol(par.buf)]))

print(Q.sd)

alpha = 1 # start flt

for (i in 1:NSTEPS) # loop over sampling first

{

old = g.BBstep(old, k, n, X, Q.sd.=Q.sd, npar, ndat) # new model with LL

# flt

ori = sample(1:nsave, 1) # the model that gets writen over

alpha = (alpha*nsave)/(alpha+nsave) # prob current model

r = runif(1) # random value

if (r <= alpha)

{

par.buf.temp = c(i+nsave, old)

par.buf[ori,] = par.buf.temp

}

}

phy = par.buf[,ncol(par.buf)]

rho = (sum(n)+phy)/(phy+1)

rownames(par.buf) = 1:nsave

#print(par.buf)

par.buf = cbind(par.buf, rho)

index.s = sort.int(par.buf[,1],index.return = TRUE)

#print(index.s)

par.buf = par.buf[index.s$ix,]

#print(par.buf)

M.list = 3:(ncol(par.buf))

SD.list = 3:(ncol(par.buf))

LB.list = 3:(ncol(par.buf))

UB.list = 3:(ncol(par.buf))

for (i in 3:(ncol(par.buf)))

{

j = i-2

M.list[j] = mean(par.buf[,i])

SD.list[j] = sd(par.buf[,i])

LB.list[j] = quantile(par.buf[,i], prob=0.025)

UB.list[j] = quantile(par.buf[,i], prob=0.975)

}

tab = data.frame(Estimate = M.list, SD = SD.list, LB = LB.list, UB = UB.list)

rownames(tab) = c(“Intercept”, colnames(X), “Psi”, “Rho”)

print(round(tab, 4))

print(round(cor(par.buf[,3:9]),4))

par(mfrow=c(1,2))

plot(par.buf[,1], phy, xlab=”Index”, ylab=”Percision”)

#quartz()

plot(par.buf[,1], rho, xlab=”Index”, ylab=”Dispersion”)

#breaks.rho = seq(1, 1.01*max(rho), length.out=201)

#h1=hist(rho, breaks=breaks.rho, plot =F)

#plot(h1$den,, h1$mids, type=”n”, xaxs=”i”, yaxs=”i”)

#polygon( c(0, h1$den, 0), c(min(h1$mids), h1$mids, max(h1$mids)), col=”grey”)

}

g.BB.glm(c(rbinom(15, size =10, prob = 0.05), rbinom(15, size =10, prob = 0.95),rbinom(30, size =10, prob = 0.5) ),

rep(10, 20),

X=as.data.frame(matrix(c(rep(0,30), rep(1,30), rnorm(60*4)), nrow=60, ncol=5)))

What does that look like? The output