inst/doc/FittingLogLinearModelsInCvam.R

In cvam: Coarsened Variable Modeling

### R code from vignette source 'FittingLogLinearModelsInCvam.Rnw'

###################################################
### code chunk number 1: "UCBA.A"
###################################################
# show the structure of the object
str(UCBAdmissions)
# display slices of the table corresponding to Dept "A" and "B"
UCBAdmissions[,,1:2]


###################################################
### code chunk number 2: "UCBA.B"
###################################################
# one-way table for Dept (dimension 3)
apply( UCBAdmissions, 3, sum )
# two-way table for Gender x Admit (dimensions 2 and 1),
# with chisquare test for independence
GenderByAdmit <- apply( UCBAdmissions, c(2,1), sum )
chisq.test( GenderByAdmit )


###################################################
### code chunk number 3: "UCBA.C"
###################################################
library(cvam)
# display the first few rows
head(microUCBAdmissions)


###################################################
### code chunk number 4: "UCBA.D"
###################################################
dF <- microUCBAdmissions  # to save typing
# this reproduces the 3-way table UCBAdmissions
result <- table( Admit = dF$Admit, 
   Gender = dF$Gender, Dept = dF$Dept )
str(result)
all.equal( result, UCBAdmissions )
# do the same thing with xtabs, which accepts formula notation
result <- xtabs( ~ Admit + Gender + Dept, data=microUCBAdmissions )


###################################################
### code chunk number 5: "UCBA.E"
###################################################
result <- as.data.frame(UCBAdmissions) 
head(result)


###################################################
### code chunk number 6: "UCBAe"
###################################################
# create a Freq variable and fill it with ones
microUCBAdmissions$Freq <- 1
# use aggregate to sum the Freq variable within categories of
# Admit, Gender, and Dept
result <- aggregate( Freq ~ Admit + Gender + Dept,
   data=microUCBAdmissions, FUN=sum )
head(result)


###################################################
### code chunk number 7: "UCBA.F"
###################################################
# from a table, specifying the row and column variables
ftable( UCBAdmissions, row.vars=c("Dept","Gender"), col.vars="Admit")
# from microdata, using a formula interface
ftable( Admit ~ Dept + Gender, data=microUCBAdmissions )
# with one row variable and two column variables 
ftable( UCBAdmissions, row.vars=c("Dept"),
   col.vars=c("Gender","Admit"))
# omitted variables are summed over
ftable( Admit ~ Gender, data=microUCBAdmissions )


###################################################
### code chunk number 8: "loglinExA"
###################################################
# display observed marginal table and odds ratio
marg <- apply( UCBAdmissions, c(2,1), sum )
marg
marg[1,1] * marg[2,2] / ( marg[2,1] * marg[1,2] )


###################################################
### code chunk number 9: "loglinExB"
###################################################
# display odds ratios for each department
UCBAdmissions[1,1,] * UCBAdmissions[2,2,] / 
   ( UCBAdmissions[1,2,] * UCBAdmissions[2,1,] )


###################################################
### code chunk number 10: "loglinExC"
###################################################
dF <- as.data.frame(UCBAdmissions)
M0 <- glm( Freq ~ Dept*Gender + Dept*Admit, family=poisson(), data=dF )
M1 <- glm( Freq ~ Dept*Gender + Dept*Admit + Gender*Admit,
   family=poisson(), data=dF )
M2 <- glm( Freq ~ Dept*Gender*Admit, family=poisson(), data=dF )


###################################################
### code chunk number 11: "loglinExD"
###################################################
dF$muHat0 <- predict(M0, type="response")
dF$muHat1 <- predict(M1, type="response")
dF$muHat2 <- predict(M2, type="response")
fit0 <- xtabs( muHat0 ~ Admit + Gender + Dept, data=dF )
fit1 <- xtabs( muHat1 ~ Admit + Gender + Dept, data=dF )
fit2 <- xtabs( muHat2 ~ Admit + Gender + Dept, data=dF )


###################################################
### code chunk number 12: "loglinExE"
###################################################
# under M0, the fitted conditional OR's should be 1.0:
fit0[1,1,] * fit0[2,2,] / ( fit0[1,2,] * fit0[2,1,] )

# under M1, the fitted conditional OR's should be equal:
fit1[1,1,] * fit1[2,2,] / ( fit1[1,2,] * fit1[2,1,] )

# under M2, the fitted conditional OR's should vary, and they
# should agree with corresponding OR's based on the observed
# frequencies, because M2 is saturated:
fit2[1,1,] * fit2[2,2,] / ( fit2[1,2,] * fit2[2,1,] )


###################################################
### code chunk number 13: "loglinExF"
###################################################
anova(M0,M1,M2)


###################################################
### code chunk number 14: "loglinExFa"
###################################################
d01 <- deviance(M0)-deviance(M1)
d12 <- deviance(M1)-deviance(M2)
d02 <- deviance(M0)-deviance(M2)


###################################################
### code chunk number 15: "loglinExGa"
###################################################
# make a list of 6 data frames, one per department
list2x2 <- as.list(1:6)
for( j in 1:6 ) list2x2[[j]] <- subset(dF, Dept==levels(dF$Dept)[j]  )
# function for computing deviance for LR test of independence
# within a department
myFunc <- function( dF ) {
   M <- glm( Freq ~ Gender + Admit, family=poisson(), data=dF )
   deviance(M)
}
# apply LR test to each department, returning a vector of deviances
dev <- sapply( list2x2, myFunc )
dev
sum(dev)


###################################################
### code chunk number 16: "loglinExH"
###################################################
# fit M0, M1, and M2 using loglin
M0 <- loglin( UCBAdmissions, margin=list( c(3,2), c(3,1) ), fit=TRUE ) 
M1 <- loglin( UCBAdmissions, margin=list( c(3,2), c(3,1), c(2,1) ), fit=TRUE ) 
M2 <- loglin( UCBAdmissions, margin=list( c(3,2,1)), fit=TRUE ) 


###################################################
### code chunk number 17: "loglinExIa"
###################################################
max( abs( fit0 - M0$fit ) )
max( abs( fit1 - M1$fit ) )
max( abs( fit2 - M2$fit ) )


###################################################
### code chunk number 18: "loglinExIb"
###################################################
M1 <- loglin( UCBAdmissions, margin=list( c(3,2), c(3,1), c(2,1) ),
   fit=TRUE, eps=1e-06 ) 
max( abs( fit1 - M1$fit ) )


###################################################
### code chunk number 19: "loglinExIc"
###################################################
M0$lrt
M1$lrt
M2$lrt


###################################################
### code chunk number 20: "loglinExJ"
###################################################
library(cvam)
dF <- as.data.frame(UCBAdmissions)
M0 <- cvam( ~ Dept*Gender + Dept*Admit, data=dF, freq=Freq )
M1 <- cvam( ~ Dept*Gender + Dept*Admit + Gender*Admit, data=dF, freq=Freq )
M2 <- cvam( ~ Dept*Gender*Admit, data=dF, freq=Freq )
anova(M0,M1,M2)


###################################################
### code chunk number 21: "loglinExK"
###################################################
get.coef(M0, withSE=TRUE)


###################################################
### code chunk number 22: "loglinExL"
###################################################
# display the fitted means for the first few cells
head( get.fitted(M0, type="mean" ) )


###################################################
### code chunk number 23: "loglinExM"
###################################################
# refit M0 with microdata to see that results are the same
M0 <- cvam( ~ Dept*Gender + Dept*Admit, data=microUCBAdmissions )
get.coef(M0, withSE=TRUE)


###################################################
### code chunk number 24: "pearsonPlotFake" (eval = FALSE)
###################################################
## fit0 <- get.fitted(M0, type="mean")
## pearson <- ( fit0$freq - fit0$fit ) / sqrt( fit0$fit )
## labs <- as.character( fit0$Dept : fit0$Gender : fit0$Admit )
## plot( 1:NROW(fit0), pearson,
##    xlab="Cell index", ylab="Pearson residual" )
## identify( 1:NROW(fit0), pearson, labels=labs )


###################################################
### code chunk number 25: FittingLogLinearModelsInCvam.Rnw:964-972
###################################################
fit0 <- get.fitted(M0, type="mean")
pearson <- ( fit0$freq - fit0$fit ) / sqrt( fit0$fit )
labs <- as.character( fit0$Dept : fit0$Gender : fit0$Admit )
plot( 1:NROW(fit0), pearson, xlab="Cell index", ylab="Pearson residual" )
w <- 7
text( w, pearson[w], labels=labs[w], pos=4 )
w <- 19
text( w, pearson[w], labels=labs[w], pos=2 )


###################################################
### code chunk number 26: "loglinExN"
###################################################
# compute the deviance for model M0
M0 <- cvam( ~ Dept*Gender + Dept*Admit, data=dF, freq=Freq )
M2 <- cvam( ~ Dept*Gender*Admit, data=dF, freq=Freq, saturated=TRUE )
dev.M0 <- -2 * ( get.loglik(M0) - get.loglik(M2) )
dev.M0


###################################################
### code chunk number 27: "loglinExO"
###################################################
# fit M1 as a conditional model
M1 <- cvam( ~ Dept*Gender + Dept*Admit + Gender*Admit | Dept + Gender,
   data=dF, freq=Freq )


###################################################
### code chunk number 28: "loglinExP"
###################################################
# show the first few fitted probabilities
head( get.fitted(M1, type="prob") )


###################################################
### code chunk number 29: "crimeA"
###################################################
data(crime)   # load the crime dataset distributed with cvam
crime
sum(crime$n)


###################################################
### code chunk number 30: "crimeB"
###################################################
# fit the model of independence
M0 <- cvam( ~ V1 + V2, freq=n, data=crime )
# fit the model of non-independence
M1 <- cvam( ~ V1 * V2, freq=n, data=crime )
# compare them
anova(M0,M1, pval=TRUE)


###################################################
### code chunk number 31: "crimeC"
###################################################
summary(M0)


###################################################
### code chunk number 32: "crimeD"
###################################################
dF <- get.fitted(M0, type="mean")
dF


###################################################
### code chunk number 33: "crimeE"
###################################################
( dF$fit[1] * dF$fit[4] ) / ( dF$fit[2] * dF$fit[3] )


###################################################
### code chunk number 34: "crimeF"
###################################################
( dF$freq[1] * dF$freq[4] ) / ( dF$freq[2] * dF$freq[3] )


###################################################
### code chunk number 35: "crimeG"
###################################################
# examine the frames from get.fitted for M0 and M1
# to make sure that they use the same cell ordering
get.fitted(M0)
get.fitted(M1)
# compute the quasi-Pearson residuals
muHat <- get.fitted(M0, type="mean")$fit
fHatSat <- get.fitted(M1, type="mean")$freq
quasiPearson <- ( fHatSat - muHat ) / sqrt( muHat )
quasiPearson


###################################################
### code chunk number 36: "raceHispA"
###################################################
data(abortion2000)
CenRace <- addNA(abortion2000$CenRace)
Hisp <- addNA(abortion2000$Hisp)
RH <- Hisp:CenRace
RH <- droplevels(RH)
levels(RH) <- list(
   nonHispWhite = "nonHisp:White",
   nonHispBlack = "nonHisp:Black",
   nonHispOther = "nonHisp:Other",
   Hisp = c("Hisp:White", "Hisp:Black", "Hisp:Hisp", "Hisp:Other", "Hisp:NA"),
   nonHispNA = "nonHisp:NA",
   NAWhite = "NA:White" )
RH  <- coarsened( RH, levelsList = list(
   nonHispNA = c("nonHispWhite", "nonHispBlack", "nonHispOther"), 
   NAWhite = c("nonHispWhite", "Hisp" ) ) )
summary(RH)


###################################################
### code chunk number 37: "raceHispB"
###################################################
baseLevels(RH)


###################################################
### code chunk number 38: "raceHispC"
###################################################
coarseLevels(RH)
mapping(RH)


###################################################
### code chunk number 39: "raceHispD"
###################################################
result <- cvam( ~ RH )
summary(result, showCoef=FALSE)
# display the fitted proportions
get.fitted(result)


###################################################
### code chunk number 40: "AbAnyA"
###################################################
# copy the four variables into a data frame
dF <- data.frame( Sex = abortion2000$Sex, RH = RH,
   PolViews = abortion2000$PolViews, AbAny = abortion2000$AbAny )
# display the first few rows
head(dF)


###################################################
### code chunk number 41: "AbAnyB"
###################################################
myFormula <- ~ Sex*RH*PolViews + AbAny*Sex + AbAny*RH + AbAny*PolViews


###################################################
### code chunk number 42: "AbAnyC"
###################################################
myMod <- cvam( myFormula, data=dF )
satMod <- cvam( ~ Sex*RH*PolViews*AbAny, data=dF, saturated=TRUE )
anova( myMod, satMod, pval=TRUE )
anova( myMod, satMod, method="AIC" )
# compute and summarize the fitted values
muHat <- get.fitted(myMod, type="mean")$fit
summary( muHat )


###################################################
### code chunk number 43: "AbAnyD"
###################################################
noSex <- cvam( ~ Sex*RH*PolViews + AbAny*RH + AbAny*PolViews, data=dF)
anova( noSex, myMod, pval=TRUE )
noRH  <- cvam( ~ Sex*RH*PolViews + AbAny*Sex + AbAny*PolViews, data=dF)
anova( noRH, myMod, pval=TRUE )
noPol <- cvam( ~ Sex*RH*PolViews + AbAny*Sex + AbAny*RH, data=dF)
anova( noPol, myMod, pval=TRUE )


###################################################
### code chunk number 44: "cvamControlA" (eval = FALSE)
###################################################
## cvamControl()


###################################################
### code chunk number 45: "cvamControlB"
###################################################
# use a boundary criterion that is less strict
satMod <- cvam( ~ Sex*RH*PolViews*AbAny, data=dF, saturated=TRUE,
   control=list(critBoundary=1e+06 ) )


###################################################
### code chunk number 46: "cvamControlC"
###################################################
round( get.fitted(satMod, type="prob")$fit, 6)


###################################################
### code chunk number 47: "priorA"
###################################################
myPrior <- cvamPrior( flatten=7.2 )


###################################################
### code chunk number 48: "priorB"
###################################################
# re-fit and compare models using the flattening constant
myMod <- cvam( myFormula, data=dF, prior=myPrior )
satMod <- cvam( ~ Sex*RH*PolViews*AbAny, data=dF,
   saturated=TRUE, prior=myPrior )


###################################################
### code chunk number 49: "priorC"
###################################################
anova( myMod, satMod, pval=TRUE, method="logP")


###################################################
### code chunk number 50: "cvamOnCvam"
###################################################
# fit the saturated model to the crime data
result <- cvam( ~ V1 * V2, data=crime, freq=n)
# run it again, starting from the previous result
result <- cvam(result)
summary(result, showCoef=FALSE)


###################################################
### code chunk number 51: "AbAnyF"
###################################################
cvamEstimate( ~ Sex + AbAny, myMod )


###################################################
### code chunk number 52: "AbAnyG"
###################################################
# estimated conditional probabilities for AbAny given Sex
cvamEstimate( ~ AbAny | Sex, myMod )


###################################################
### code chunk number 53: "AbAnyH"
###################################################
# conditional probabilities for AbAny given RH and PolViews
est <- cvamEstimate( ~ AbAny | RH + PolViews, myMod ) 

# reshape the probabilities into a three-dimensional array
xtab <- xtabs( prob ~ AbAny + RH + PolViews, data = est )

# display the array as a flat table
ftable( xtab, row.vars=c("PolViews", "RH"), col.vars="AbAny" )


###################################################
### code chunk number 54: "predictA"
###################################################
# display the crime data
crime
# fit the model of non-independence
fit <- cvam( ~ V1 * V2, data=crime, freq=n )
# display predictions for V1
cvamPredict( ~ V1, fit, data=crime )


###################################################
### code chunk number 55: "predictB"
###################################################
# display predicted frequencies for V1
cvamPredict( ~ V1, fit, data=crime, freq=n )


###################################################
### code chunk number 56: "predictC"
###################################################
# display predicted frequencies for V1 and V2
cvamPredict( ~ V1 + V2, fit, data=crime, freq=n )


###################################################
### code chunk number 57: "imputeA"
###################################################
set.seed(69852)
cvamImpute( fit, data=crime )


###################################################
### code chunk number 58: "imputeB"
###################################################
cvamImpute( fit, data=crime, freq=n )


###################################################
### code chunk number 59: "approxBayesA"
###################################################
# fit the non-independence model to the crime data
fitML <- cvam( ~ V1 * V2, data=crime, freq=n )
# display the ML estimate for beta and pi
get.coef( fitML )
get.fitted( fitML, type="prob" )$fit
# draw from the approximate posterior, display new beta and pi
set.seed(83425)
obj <- cvam(fitML, method="approxBayes")
get.coef( obj )
get.fitted( obj, type="prob" )$fit


###################################################
### code chunk number 60: "approxBayesB"
###################################################
# produce 5,000 draws of beta, saving also the resulting pi vectors
obj <- cvam(fitML, method="approxBayes",
   control=list(iterApproxBayes=5000, saveProbSeries=TRUE) )
# display the first few beta and pi vectors
head( get.coefSeries(obj) )
head( get.probSeries(obj) )


###################################################
### code chunk number 61: "approxBayesC"
###################################################
pi.series <- get.probSeries(obj)
delta <- pi.series[,3] - pi.series[,2]
summary(delta)
sum( delta > 0 )


###################################################
### code chunk number 62: "mcmcA"
###################################################
set.seed(4358)
fit <- cvam( ~ V1 * V2, data=crime, freq=n, method="MCMC")
summary(fit)


###################################################
### code chunk number 63: "mcmcB"
###################################################
betaSeries <- get.coefSeries( fit )
library(coda)
summary( betaSeries )


###################################################
### code chunk number 64: "mcmcC"
###################################################
# display trace plots and density estimates
plot( betaSeries )


###################################################
### code chunk number 65: FittingLogLinearModelsInCvam.Rnw:1896-1897
###################################################
plot( betaSeries )


###################################################
### code chunk number 66: "mcmcE"
###################################################
# display autocorrelation plots
acfplot( betaSeries )


###################################################
### code chunk number 67: FittingLogLinearModelsInCvam.Rnw:1911-1912
###################################################
acfplot( betaSeries )


###################################################
### code chunk number 68: "mcmcG"
###################################################
get.fitted(fit, type="prob")


###################################################
### code chunk number 69: "mcmcg"
###################################################
set.seed(4358)
fit <- cvam( ~ V1 * V2, data=crime, freq=n, method="MCMC",
   control=list( saveProbSeries=TRUE ) )
piSeries <- get.probSeries(fit)
delta <- piSeries[,3] - piSeries[,2]
summary(delta)
sum( delta > 0 )


###################################################
### code chunk number 70: "miA"
###################################################
impList <- as.list(1:10) # a list to store the imputed datasets
set.seed(769090)         # for reproducibility
for(m in 1:10) {
   # run MCMC under the non-independence model
   tmp <- cvam( ~ V1 * V2, data=crime, freq=n, method="MCMC")
   # impute under the simulated parameters
   impList[[m]] <- cvamImpute( tmp, crime, freq=n)
}
# display the first two imputations
impList[1:2]


###################################################
### code chunk number 71: "miB"
###################################################
# run MCMC for 5,000 iterations, saving an imputation at every 500th
result <- cvam( ~ V1 * V2, data=crime, freq=n, method="MCMC",
   control=list( iterMCMC=5000, imputeEvery=500 ) )
get.imputedFreq(result)


###################################################
### code chunk number 72: "miC"
###################################################
#  run EM, then create ten imputations with approxBayes
fitML <- cvam( ~ V1 * V2, data=crime, freq=n ) 
result <- cvam( fitML, method="approxBayes",
   control=list( iterApproxBayes=10, imputeApproxBayes=TRUE ) )
get.imputedFreq(result)


###################################################
### code chunk number 73: "miD"
###################################################
set.seed(54981)
result <- cvam( fitML, method="MCMC",
   control=list( iterMCMC=5000, imputeEvery=500 ) )
impData <- get.imputedFreq(result)[-(1:2)] # just the frequencies 
est.list <- std.err.list <- as.list(1:10)  # to hold the estimates and SEs
for( m in 1:10 ) {
   f <- impData[,m]
   est.list[[m]] <- log( (f[1] * f[4]) / (f[2] * f[3]) )
   std.err.list[[m]] <- sqrt( sum(1/f) )
}
miInference( est.list, std.err.list )


###################################################
### code chunk number 74: "micromiA"
###################################################
# put the four variables into a data frame
dF <- data.frame( Sex = abortion2000$Sex, RH = RH,
   PolViews = abortion2000$PolViews, AbAny = abortion2000$AbAny )
# fit the saturated model with EM, then do a test run of MCMC
fitEM <- cvam( ~ Sex * RH * PolViews * AbAny, data=dF )
set.seed(598902)
fitMCMC <- cvam( fitEM, method="MCMC")


###################################################
### code chunk number 75: "micromiB"
###################################################
# display fitted cell probs, rounded to five decimal places
round( get.fitted(fitEM, type="prob", mfTrue=FALSE ), 5)
# display some of the coefs and SEs
head( get.coef(fitEM, withSE=TRUE) )


###################################################
### code chunk number 76: "micromiC"
###################################################
# re-run EM with a ridge factor of 0.5
fitEM.ridge <- cvam( ~ Sex * RH * PolViews * AbAny, data=dF,
   prior=cvamPrior( ridge=.5 ) )
round( get.fitted(fitEM.ridge, type="prob", mfTrue=FALSE ), 5)
head( get.coef(fitEM.ridge, withSE=TRUE) )


###################################################
### code chunk number 77: "micromiD"
###################################################
-2 * ( get.loglik(fitEM.ridge) - get.loglik(fitEM) )
exp( get.loglik(fitEM) - get.loglik(fitEM.ridge) )


###################################################
### code chunk number 78: "micromiE"
###################################################
set.seed(87900)
fitMCMC <- cvam( fitEM.ridge, method="MCMC" )


###################################################
### code chunk number 79: "micromiF"
###################################################
set.seed(87900)
fitMCMC <- cvam( fitEM.ridge, method="MCMC",
   control=list( typeMCMC="RWM", tuneRWM=c(1000,.17) ) )


###################################################
### code chunk number 80: "micromiG"
###################################################
M <- 25
impList <- as.list(1:M)  # dummy list to hold the imputed datasets
set.seed(2343)
for( m in 1:M ) {
   # take 2,500 steps of MCMC, then impute
   fitMCMC <- cvam( fitMCMC, control=list(iterMCMC=2500) )
   impList[[m]] <- cvamImpute( fitMCMC, data=dF )
}
# display the first few rows of the original data and 
# the first imputed dataset
head( dF )
head( impList[[1]] )


###################################################
### code chunk number 81: "micromiH"
###################################################
est.list <- SE.list <- as.list(1:M)
for( m in 1:M ) {
   # extract the imputed dataset
   impData <- impList[[m]]
   # create the binary response and fit the logit model
   impData$y <- 1 * ( impData$AbAny == "Yes" )
   logitFit <- glm( y ~ Sex + RH + PolViews, data=impData,
      family=binomial() )
   # extract matrix of coefficients and SEs
   coefMat <- summary(logitFit)$coef
   est.list[[m]] <- coefMat[,1]
   SE.list[[m]]  <- coefMat[,2]
}
# combine the results with Rubin's rules
miInference( est.list, SE.list )


###################################################
### code chunk number 82: "syntheticA"
###################################################
# take 2,500 more steps of MCMC and draw a synthetic dataset
fitMCMC <- cvam( fitMCMC )
synthData <- cvamImpute( fitMCMC, data=dF)
head( synthData )


###################################################
### code chunk number 83: "hivA"
###################################################
hivtest


###################################################
### code chunk number 84: "hivB"
###################################################
hivtest$L <- latentFactor( NROW(hivtest), 2 )
hivtest


###################################################
### code chunk number 85: "hivC"
###################################################
# set the RNG seed and fit the model of local independence
set.seed(125)
fit <- cvam( ~ L*A + L*B + L*C + L*D, data=hivtest, freq=COUNT,
   control = list( startValJitter=.1 ) )


###################################################
### code chunk number 86: "hivE"
###################################################
cvamEstimate( list( ~L, ~A|L, ~B|L, ~C|L, ~D|L ), fit )


###################################################
### code chunk number 87: "hivFa"
###################################################
# perform the lack-of-fit test
fitSat <- cvam( ~ A*B*C*D, data=hivtest, freq=COUNT )
anova( fit, fitSat, pval=TRUE )


###################################################
### code chunk number 88: "hivFb"
###################################################
satFrame <- get.fitted( fitSat, type="mean" )
# this frame has 16 rows; display the first few
head(satFrame)
# get rid of the fitted values, because they are redundant
satFrame$fit <- NULL


###################################################
### code chunk number 89: "hivFbb"
###################################################
LCFrame <-  get.fitted( fit, type="mean" )
# this frame has 32 rows; display the first few
head(LCFrame)


###################################################
### code chunk number 90: "hivFc"
###################################################
muHatTable <- xtabs( fit ~ A + B + C + D, data=LCFrame )
muHatFrame <- as.data.frame( muHatTable, responseName = "muHat" )
# diplay the first few rows to make sure that the
# cell order is correct
head( muHatFrame)


###################################################
### code chunk number 91: "hivFd"
###################################################
muHat <- muHatFrame$muHat
quasiPearson <- ( satFrame$freq - muHat ) / sqrt( muHat )
satFrame$muHat <- round( muHat, 3 )
satFrame$quasiPearson <- round( quasiPearson, 2 )
satFrame


###################################################
### code chunk number 92: "hivI"
###################################################
set.seed(85657)
fitLAB <- cvam( ~ L*A + L*B + L*C + L*D + L*A*B, 
   data=hivtest, freq=COUNT,
   control = list(startValJitter=.1) )
anova(fit, fitLAB, fitSat, pval=TRUE)
fitLAC <- cvam( ~ L*A + L*B + L*C + L*D + L*A*C, 
   data=hivtest, freq=COUNT,
   control = list(startValJitter=.1) )
anova(fit, fitLAC, fitSat, pval=TRUE)
fitLAD <- cvam( ~ L*A + L*B + L*C + L*D + L*A*D, 
   data=hivtest, freq=COUNT,
   control = list(startValJitter=.1) )
anova(fit, fitLAD, fitSat, pval=TRUE)
fitLBC <- cvam( ~ L*A + L*B + L*C + L*D + L*B*C, 
   data=hivtest, freq=COUNT,
   control = list(startValJitter=.1) )
anova(fit, fitLBC, fitSat, pval=TRUE)
fitLBD <- cvam( ~ L*A + L*B + L*C + L*D + L*B*D, 
   data=hivtest, freq=COUNT,
   control = list(startValJitter=.1) )
anova(fit, fitLBD, fitSat, pval=TRUE)
fitLCD <- cvam( ~ L*A + L*B + L*C + L*D + L*C*D, 
   data=hivtest, freq=COUNT,
   control = list(startValJitter=.1) )
anova(fit, fitLCD, fitSat, pval=TRUE)


###################################################
### code chunk number 93: "hivJ"
###################################################
fitBoth <- cvam( ~ L*A + L*B + L*C + L*D + L*A*D + L*B*C, 
   data=hivtest, freq=COUNT,
   control = list(startValJitter=.1) )
anova(fitLAD, fitBoth)
anova(fitLBC, fitBoth)


###################################################
### code chunk number 94: hivK
###################################################
# get predicted probabilities and display them with the dataset
pred <- cvamPredict( ~L, fitLBC, data=hivtest )
cbind( hivtest, round(pred, 3) )


###################################################
### code chunk number 95: "hivk"
###################################################
predFrame <- hivtest[1:8,]
predFrame$COUNT <- NULL
predFrame[["A"]][] <- NA
predFrame[["B"]][] <- NA
predFrame[["C"]][] <- NA
predFrame[["D"]][] <- NA
predFrame[["A"]][1] <- "pos"; predFrame[["A"]][2] <- "neg"
predFrame[["B"]][3] <- "pos"; predFrame[["B"]][4] <- "neg"
predFrame[["C"]][5] <- "pos"; predFrame[["C"]][6] <- "neg"
predFrame[["D"]][7] <- "pos"; predFrame[["D"]][8] <- "neg"
predFrame[["A"]] <- coarsened( predFrame[["A"]] )
predFrame[["B"]] <- coarsened( predFrame[["B"]] )
predFrame[["C"]] <- coarsened( predFrame[["C"]] )
predFrame[["D"]] <- coarsened( predFrame[["D"]] )
pred <- cvamPredict( ~L, fitLBC, data=predFrame )
cbind( predFrame, round(pred, 3) )


###################################################
### code chunk number 96: "hivL"
###################################################
pred <- cvamPredict( ~L, fit, data=predFrame )
cbind( predFrame, round(pred, 3) )


###################################################
### code chunk number 97: "hivM"
###################################################
# re-fit the model with EM using a small ridge factor
set.seed(7666)
fitLBC <- cvam( ~ L*A + L*B + L*C + L*D + L*B*C, 
   data=hivtest, freq=COUNT, prior=cvamPrior( ridge=.1 ),
   control = list(startValJitter=.1) )


###################################################
### code chunk number 98: "hivMa"
###################################################
# do a long run of MCMC and save ten imputed datasets
fitMCMC <- cvam(fitLBC, method="MCMC",
   control=list( typeMCMC="RWM", tuneRWM=c(1000,.5),
      iterMCMC=25000, imputeEvery=2500 ) )


###################################################
### code chunk number 99: "hivN"
###################################################
# check to see if any label switching has occurred
impData <- get.imputedFreq(fitMCMC)
head(impData)


###################################################
### code chunk number 100: "hivNa"
###################################################
impData$freq <- impData[["imp.1"]] # first imputation
BCL <- xtabs( freq ~ B + C + L, data=impData )
BCL


###################################################
### code chunk number 101: "hivO"
###################################################
# use multiple imputations to examine the conditional
# BC odds ratios given L=1 and L=2
est.list <- SE.list <- as.list(1:10)
for( m in 1:10 ) {
   # get the imputed marginal table BxCxL
   impName <- paste( "imp", format(m), sep="." )
   impData$freq <- impData[[impName]]
   BCL <- xtabs( freq ~ B + C + L, data=impData )
   # add 1/2 to every cell to avoid problems
   BCL <- BCL + .5
   # get BC log-odds ratio and SE for L=1
   BCL.1 <- BCL[,,"1"]
   logOR.1 <- log( ( BCL.1[1,1] * BCL.1[2,2] ) /
      ( BCL.1[1,2] * BCL.1[2,1] ) )
   SE.1 <- sqrt( sum( 1/BCL.1 ) )
   # get BC log-odds ratio and SE for L=2
   BCL.2 <- BCL[,,"2"]
   logOR.2 <- log( ( BCL.2[1,1] * BCL.2[2,2] ) /
      ( BCL.2[1,2] * BCL.2[2,1] ) )
   SE.2 <- sqrt( sum( 1/BCL.2 ) )
   # save the estimates and SEs
   est.list[[m]] <- c( logOR.1, logOR.2 )
   SE.list[[m]] <- c( SE.1, SE.2 )
}
miInference( est.list, SE.list )