# ----------------------------------------------------------------------
#
# PURPOSE:  Compute the relative efficiency for overdispersed count data
#
# AUTHOR:  P. Heagerty/D. Gillen
#
# DATE:  27Jan2002
#        02Feb2007
#
# ----------------------------------------------------------------------
#
#####  Function to sample data
#
getData <- function( beta ){
   nobs <- 200
   #
   ##### a linear predictor
   #
   x1 <- (c(1:nobs)-nobs/2)/(nobs/2)
   #
   ##### a binary predictor
   #
   x2 <- rep( c(0,1), c(nobs,nobs)/2 )
   #
   X <- cbind( 1, x1, x2 )
   mu <- exp( X %*% beta )
   data <- list( X, mu )
   names( data ) <- c( "X", "mu" )
   data
}
#
#####
#####  Part (a):  ARE when true model is N-B
#####
#
#####  Function to calculate ARE when true model is N-B
#
areNB <- function( beta, alpha ){
   data <- getData( beta )
   X <- as.matrix( data$X )
   mu <- as.vector( data$mu )
   A.poisson <- t(X)%*%diag(mu)%*%X
   B.true <- t(X)%*%diag( mu*(1+alpha*mu) )%*%X
   A.negbin <- t(X)%*%diag( mu/(1+alpha*mu) )%*%X
   #
   se.poisson <- sqrt( diag( solve(A.poisson)%*%B.true%*%solve(A.poisson) ) )
   se.optimal <- sqrt( diag( solve(A.negbin) ) )
   #
   cbind( se.poisson, se.optimal, (se.optimal)^2 / (se.poisson)^2 )
}
#
areNB( beta=c( 1, 1, .5 ), alpha=.25 ) 
areNB( beta=c( 1, 1, .5 ), alpha=.50 ) 
areNB( beta=c( 1, 1, .5 ), alpha=1.0 ) 
#
##### 
##### Part (b):  Effect of changing beta
##### 
#
areNB( beta=c( 1, 0.0, .5 ), alpha=1.0 ) 
areNB( beta=c( 1, 0.5, .5 ), alpha=1.0 ) 
areNB( beta=c( 1, 1.0, .5 ), alpha=1.0 ) 
#
#####
##### Part (c): What if really scale and we fit negbin? 
#####
#
#####  Function to calculate ARE when true model is scale, but N-B is fit
#
areScale <- function( beta, scale ){
   data <- getData( beta )
   X <- as.matrix( data$X )
   mu <- as.vector( data$mu )
   A.poisson <- t(X)%*%diag(mu)%*%X
   #
   alpha <- mean( (scale-1)/mu )
   alpha <- (scale-1)/mean( mu )
   A.negbin <- t(X)%*%diag( mu/(1+alpha*mu) )%*%X
   B.true <- t(X)%*%diag( scale*mu/(1+alpha*mu)^2 )%*%X
   se.negbin <- sqrt( diag( solve(A.negbin)%*%B.true%*%solve(A.negbin) ) )
   se.optimal <- sqrt( diag( scale* solve(A.poisson) ) )
   #
   cbind( se.negbin, se.optimal, (se.optimal)^2 / (se.negbin)^2 )
}
#
areScale( beta=c( 1, 1, .5 ), scale=2.0 ) 
areScale( beta=c( 1, 1, .5 ), scale=3.0 ) 
areScale( beta=c( 1, 1, .5 ), scale=4.0 ) 
#
##### 
##### Part (d):  Verification via simulation
##### 
#
####### (1) alpha = 1.0
#
nobs <- 200
#
##### a linear predictor
#
x1 <- (c(1:nobs)-nobs/2)/(nobs/2)
#
##### a binary predictor
#
x2 <- rep( c(0,1), c(nobs,nobs)/2 )
#
mu <- exp( 1.0 + 1.0*x1 + 0.5*x2 )
alpha <- 1.0
#
Nsim <- 1000
beta0 <- matrix( NA, Nsim, 3 )
beta1 <- matrix( NA, Nsim, 3 )
beta2 <- matrix( NA, Nsim, 3 )
#
library(MASS)
for( j in 1:Nsim ){
#
##### generate data
#
bbb <- rgamma( rep(1,nobs), shape=1/alpha ) * alpha
y <- rpois( rep(1,nobs), lambda=mu*bbb )
#
fit0 <- glm( y ~ x1 + x2, family=poisson )
beta0[j,] <- fit0$coef
#
####
#### Using GLM and weights
####
#
maxits <- 10
converge <- F
fit <- fit0
count <- 0
while( !converge & count<maxits ){
  old.beta <- fit$coef
  mmm <- fitted(fit)
  r2 <- (y-mmm)^2 / mmm
  www <- 1/(1+alpha*mmm)
  fit <- glm( y ~ x1 + x2, family=poisson, weight=www )
  new.beta <- fit$coef
  if( max( abs( new.beta-old.beta) ) < 1e-5 ) converge <- T
  count <- count + 1 
 }
beta1[j,] <- new.beta
#
#####
##### Using NB2 MLE
#####
#
simdata <- data.frame( y=y, x1=x1, x2=x2 )
fit2 <- glm.nb( y ~ x1 + x2, data=simdata )
beta2[j,] <- fit2$coef
#
}
#
s0 <- sqrt( diag( var( beta0 ) ) ) 
s1 <- sqrt( diag( var( beta1 ) ) )
s2 <- sqrt( diag( var( beta2 ) ) )
print( cbind( s0, s1, s2, (s1/s0)^2, (s2/s0)^2 ) )
#
####### (2) scale = 4.0
#
scale <- 4.0
#
Nsim <- 1000
beta0 <- matrix( NA, Nsim, 3 )
beta1 <- matrix( NA, Nsim, 3 )
beta2 <- matrix( NA, Nsim, 3 )
#
for( j in 1:Nsim ){
#
##### generate data
#
bbb <- rgamma( rep(1,nobs), shape=mu/(scale-1) ) * ( (scale-1)/mu )
y <- rpois( rep(1,nobs), lambda=mu*bbb )
#
fit0 <- glm( y ~ x1 + x2, family=poisson )
mu0 <- fitted(fit0)
beta0[j,] <- fit0$coef
#
####
#### Using GLM and weights
####
#
maxits <- 10
converge <- F
fit <- fit0
count <- 0
while( !converge & count<maxits ){
  old.beta <- fit$coef
  mmm <- fitted(fit)
  aaa <- mean( (scale - 1)/mu0 )
  www <- 1/(1+aaa*mmm)
  fit <- glm( y ~ x1 + x2, family=poisson, weight=www )
  new.beta <- fit$coef
  if( max( abs( new.beta-old.beta) ) < 1e-5 ) converge <- T
  count <- count + 1 
 }
beta1[j,] <- new.beta
#
#####
##### Using NB2 MLE
#####
#
simdata <- data.frame( y=y, x1=x1, x2=x2 )
fit2 <- glm.nb( y ~ x1 + x2, data=simdata )
beta2[j,] <- fit2$coef
#
}
#
s0 <- sqrt( diag( var( beta0 ) ) ) 
s1 <- sqrt( diag( var( beta1 ) ) )
s2 <- sqrt( diag( var( beta2 ) ) )
print( cbind( s0, s1, s2, (s0/s1)^2, (s0/s2)^2 ) )
#
#####
##### Part (e):  Effect of estimation of alpha
#####
#
alpha <- 1.0
#
Nsim <- 1000
beta0 <- matrix( NA, Nsim, 3 )
scale0 <- rep( NA, Nsim )
beta1 <- matrix( NA, Nsim, 3 )
alpha1 <- rep( NA, Nsim )
beta2 <- matrix( NA, Nsim, 3 )
alpha2 <- rep( NA, Nsim )
#
for( j in 1:Nsim ){
#
##### generate data
#
bbb <- rgamma( rep(1,nobs), shape=1/alpha ) * alpha
y <- rpois( rep(1,nobs), lambda=mu*bbb )
#
fit0 <- glm( y ~ x1 + x2, family=poisson )
mu0 <- fitted(fit0)
beta0[j,] <- fit0$coef
r2 <- (y-mu0)^2 / mu0 
scale0[j] <- mean(r2)
#
####
#### Using GLM and weights
####
#
maxits <- 10
converge <- F
aaa <- mean( (r2 - 1)/mu0 )
fit <- fit0
count <- 0
while( !converge & count<maxits ){
  old.beta <- fit$coef
  mmm <- fitted(fit)
  r2 <- (y-mmm)^2 / mmm
  aaa <- mean( (r2 - 1)/mu0 )
  www <- 1/(1+aaa*mmm)
  fit <- glm( y ~ x1 + x2, family=poisson, weight=www )
  new.beta <- fit$coef
  if( max( abs( new.beta-old.beta) ) < 1e-5 ) converge <- T
  count <- count + 1 
 }
beta1[j,] <- new.beta
alpha1[j] <- aaa
#
#####
##### Using NB2 MLE
#####
#
simdata <- data.frame( y=y, x1=x1, x2=x2 )
fit2 <- glm.nb( y ~ x1 + x2, data=simdata )
beta2[j,] <- fit2$coef
alpha2[j] <- 1/fit2$theta
}
#
s0 <- sqrt( diag( var( beta0 ) ) ) 
s1 <- sqrt( diag( var( beta1 ) ) )
s2 <- sqrt( diag( var( beta2 ) ) )
print( cbind( s0, s1, s2, (s1/s0)^2, (s2/s0)^2 ) )
#
#
# end-of-file...