demo/Heijden2003.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/demo/Heijden2003.R
##
## Reproducing results from data anlysis in the study of Heijden 2003.

library(ipeglim)
rm(list=ls())
data(IINEE)
cat("This reproduction is based on the study of Heijden (2003)\n")
cat("http://stat.uibk.ac.at/smij/3.4-VanDerHeijdenEtAl-Data_and_Software.zip\n")

iinee <- IINEE

## check the data structure
head(iinee)

## table 2 (p.309)
iinee$nation0 <- with(iinee, ifelse(nation1==0 & nation2==0 & nation3==0 & nation4==0 & nation5==0, 1, 0))
tb <- with(iinee, list( 
  xtabs(~age+capture), 
  xtabs(~gender+capture), 
  xtabs(~nation1+capture)[1,], 
  xtabs(~nation2+capture)[1,], 
  xtabs(~nation3+capture)[1,],
  xtabs(~nation4+capture)[1,], 
  xtabs(~nation5+capture)[1,],
  xtabs(~nation0+capture)[2,], 
  xtabs(~reason+capture)))
tb2 <- as.table(do.call(rbind, tb))
colnames(tb2) <- paste("f", 1:6, sep="")
rownames(tb2) <- c(">40 years", "<40 years", "Female", "Male", "Turkey", "North Africa", "Rest of Africa", "Surinam", "Asia", "America, America", "Being illegal", "Other reason")

## table 3. simulation study (in another file)

## table 4. (p.318)
fit <- ztpr(formula=capture~gender+age+nation1+nation2+nation3+nation4+nation5+reason, data=iinee, dist="poisson", ztrunc=TRUE)
tb4 <- summary(fit, HT.est=TRUE, LM.test=TRUE)
tb4

## table 5. (p.319)
fit0 <- summary(ztpr(formula=capture~1, data=iinee, dist="poisson", ztrunc=TRUE), HT.est=TRUE, LM.test=TRUE)
fit1 <- summary(ztpr(formula=capture~gender, data=iinee, dist="poisson", ztrunc=TRUE), HT.est=TRUE, LM.test=TRUE)
fit2 <- summary(ztpr(formula=capture~gender+age, data=iinee, dist="poisson", ztrunc=TRUE), HT.est=TRUE, LM.test=TRUE)
fit3 <- summary(ztpr(formula=capture~gender+age+nation1+nation2+nation3+nation4+nation5, data=iinee, dist="poisson", ztrunc=TRUE), HT.est=TRUE, LM.test=TRUE)
fit4 <- summary(ztpr(formula=capture~gender+age+nation1+nation2+nation3+nation4+nation5+reason, data=iinee, dist="poisson", ztrunc=TRUE), HT.est=TRUE, LM.test=TRUE)
fit <- list(fit0, fit1, fit2, fit3, fit4)

aic <- do.call(c, lapply(fit, with, aic))
k <- do.call(c, lapply(fit, with, df)) # aic = 2*k - 2llk
G2 <- c(NA, -(diff(aic)-2*diff(k)))
df <- c(NA, diff(k))
pstar <- pchisq(q=G2, df=df, lower.tail=FALSE)
chi2 <- do.call(c, lapply(fit, with, LM.chisq))
Nhat <- do.call(c, lapply(fit, with, N))
cil <- do.call(c, lapply(fit, with, cil))
ciu <- do.call(c, lapply(fit, with, ciu))
tb5 <- cbind(aic, G2, df, pstar, chi2, Nhat, cil, ciu)
colnames(tb5) <- c("AIC", "G2", "df", "P*", "chi2", "N", "CIL", "CIU")
rownames(tb5) <- c("Null", "G", "G+A", "G+A+N", "G+A+N+R")
round(tb5,3)


# table 6. (p.319)
mu <- fit3$mu
freq <- table(iinee$capture)
k <- as.numeric(names(freq))
expected <- numeric(length(k))
for(i in k) expected[i] <- sum(dztpois(i, lambda=mu))
obs <- c(0, freq)
est <- c(fit3$N-sum(expected), expected)
res <- c(NA, (freq-expected)/sqrt(expected))
tb6 <- cbind(obs, est, res)
colnames(tb6) <- c("Observed", "Estimated", "Residuals")
rownames(tb6) <- seq(0,6,1)
tb6

# table 7. (p.320)
m3 <- ztpr(formula=capture~gender+age+nation1+nation2+nation3+nation4+nation5, data=iinee, dist="poisson", ztrunc=TRUE)
X <- model.matrix(obj=capture~gender+age+nation1+nation2+nation3+nation4+nation5, data=iinee)

m3$X <- X[which(iinee$gender==1),]
sub3.11 <- summary(m3, HT.est=TRUE, LM.test=FALSE)
m3$X <- X[which(iinee$gender==0),]
sub3.10 <- summary(m3, HT.est=TRUE, LM.test=FALSE)

m3$X <- X[which(iinee$age==1),]
sub3.21 <- summary(m3, HT.est=TRUE, LM.test=FALSE)
m3$X <- X[which(iinee$age==0),]
sub3.20 <- summary(m3, HT.est=TRUE, LM.test=FALSE)

m3$X <- X[which(iinee$nation1==1),]
sub3.31 <- summary(m3, HT.est=TRUE, LM.test=FALSE)

m3$X <- X[which(iinee$nation2==1),]
sub3.41 <- summary(m3, HT.est=TRUE, LM.test=FALSE)

m3$X <- X[which(iinee$nation3==1),]
sub3.51 <- summary(m3, HT.est=TRUE, LM.test=FALSE)

m3$X <- X[which(iinee$nation4==1),]
sub3.61 <- summary(m3, HT.est=TRUE, LM.test=FALSE)

m3$X <- X[which(iinee$nation5==1),]
sub3.71 <- summary(m3, HT.est=TRUE, LM.test=FALSE)

sub <- list(sub3.11, sub3.10, sub3.21, sub3.20, sub3.31, sub3.41, 
  sub3.51, sub3.61, sub3.71)

n <- do.call(rbind, lapply(sub, with, n))
N <- do.call(rbind, lapply(sub, with, N))
cil <- do.call(rbind, lapply(sub, with, cil))
ciu <- do.call(rbind, lapply(sub, with, ciu))

tb7 <- cbind(n, N, cil, ciu, n/N)
colnames(tb7) <- c("Observed", "Expected", "cil", "ciu", "Rate")
rownames(tb7) <- c("Male", "Female", "Age<40", "Age>40", "Turkey", "N.Africa", "R.Africa", "Surinam", "Asia")
tb7