#' Fit a hierarchical log-linear model with algebraic methods
#'
#' [loglinear()] fits a hierarchical log-linear model to a dataset (typically a
#' contingency table) and performs an exact conditional test on the fitted model
#' using various distance metrics. The exact test, which is a goodness-of-fit
#' test, is performed via Monte-Carlo sampling from the conditional distribution
#' of the table given the sufficient statistics of the model. In short,
#' inference is drawn by comparing the statistic of the observed table (be it an
#' unnormalized log-likelihood, a Pearson chi-squared value, or another metric)
#' to those of the samples. The proportion of sampled tables with equal to or
#' more extreme values than the observed table is the resulting p-value.
#'
#' In many ways, [loglinear()] is like [stats::loglin()] or [MASS::loglm()];
#' however, there are a few key differences.
#'
#' The first difference is methodological. The tests conducted with
#' [loglinear()] are exact conditional tests based on the conditional
#' distribution of the data given the sufficient statistics for the model. In
#' other words, they are analogues of Fisher's exact test for generic log-linear
#' models. These tests are made possible by advances in algebraic statistics;
#' see references 1--3 below.
#'
#' The second difference between [loglinear()] and [loglin()] or [loglm()] is
#' that inference is made through Monte Carlo simulation. In particular,
#' [loglinear()] leverages Markov moves to sample from the conditional
#' distribution of the data given its sufficent statistics. If the software
#' 4ti2 is installed on your machine, you can use [markov()] (or let
#' [loglinear()] use [markov()]) to generate a Markov basis to use for the
#' Markov moves. This basis is guaranteed to produce a MCMC routine that
#' converges to the conditional distribution of interest. Since [loglinear()]
#' uses Monte Carlo simulation to conduct inference, and since it uses MCMC to
#' do so, concerns typical to MCMC should be addressed. In particular, issues
#' such as burn in and mixing (autocorrelation of samples) should be addressed.
#' The examples illustrate some of these topics. The result is a p-value that
#' is generated by Monte Carlo simulation. Its standard error is provided
#' (computed as in the standard CLT confidence interval) to give a sense of the
#' Monte Carlo error.
#'
#' A third way that [loglinear()] differs from [loglin()] or [loglm()] is in
#' generalizing the kinds of tests performed. While those allow for asymptotic
#' unconditional testing using Pearson's X^2 test and the likelihood ratio test,
#' [loglinear()] gives several test statistics: Pearson's X^2, the likelihood
#' ratio G^2, Freeman-Tukey, Cressie-Read (lambda = 2/3), and Neyman's modified
#' X^2., see the last reference. In other words, to compute the exact p-value,
#' iter = 1e4 samples are sampled from the conditional distribution of the table
#' given the sufficient statistics, and then the proportion of tables that have
#' X^2, G^2, etc. values greater than or equal to that of the observed table is
#' the p value for the (conditional) exact test. A similar, and perhaps
#' preferable approach, simply adds up the probabilities of the tables that have
#' probabilities less than or equal to that of the observed table; this is the
#' first line output in hierarchical and does not use a test statistic.
#'
#' Some authors (see the third reference) suggest that for discrete problems, a
#' "mid p value" is preferable to the traditional p value, and when presented
#' should be interepreted in the same way. If the p value is defined to be
#' P(samps >= obs), the mid p value is defined to be P(samps > obs) + P(samps ==
#' obs)/2. The mid p value is computed for each test.
#'
#' @param model hierarchical log-linear model specification
#' @param data data, typically as a table but can be in different formats. see
#' [teshape()]
#' @param init the initialization of the chain. by default, this is the observed
#' table
#' @param iter number of chain iterations
#' @param burn burn-in
#' @param thin thinning
#' @param engine \code{"C++"} or \code{"R"}? (C++ is significantly faster)
#' @param method should the expected value (exp) be fit using iterative
#' proportional fitting (via loglin) or the MCMC as the average of the steps?
#' @param moves the markov moves for the mcmc (as columns of a matrix).
#' @param ... ...
#' @return a list containing named elements
#'
#' \itemize{
#'
#' \item \code{steps}: an integer matrix whose columns represent individual
#' samples from the mcmc.
#'
#' \item \code{moves}: the moves used for the proposal distribution in the
#' mcmc, computed with 4ti2 (note that only the positive moves are given).
#'
#' \item \code{accept_prob}: the empirical transition probability of the
#' moves, including the thinned moves.
#'
#' \item \code{param}: the fitted parameters of the log linear model.
#'
#' \item \code{df}: parameters per term in the model
#'
#' \item \code{quality}: model selection statistics AIC, AICc, and BIC.
#'
#' \item \code{residuals}: the (unstandardized) pearson residuals (O - E) /
#' sqrt(E)
#'
#' \item \code{call}: the call.
#'
#' \item \code{obs}: the contingency table given.
#'
#' \item \code{exp}: the fit contingency table as an integer array.
#'
#' \item \code{A}: the sufficient statistics computing matrix (from Tmaker).
#'
#' \item \code{p.value}: the exact p-values of individual tests, accurate to
#' Monte-Carlo error. these are computed as the proportion of samples with
#' statistics equal to or larger than the oberved statistic.
#'
#' \item \code{mid.p.value}: the mid p.values, see Agresti pp.20--21.
#'
#' \item \code{statistic}: the pearson's chi-squared (X2), likelihood ratio
#' (G2), Freeman-Tukey (FT), Cressie-Read (CR), and Neyman modified
#' chi-squared (NM) statistics computed for the table given.
#'
#' \item \code{sampsStats}: the statistics computed for each mcmc sample.
#'
#' \item \code{cells}: the number of cells in the table.
#'
#' \item \code{method}: the method used to estimate the table.
#'
#' }
#' @export loglinear
#' @aliases hierarchical
#' @author David Kahle
#' @seealso [stats::loglin()], [MASS::loglm()], [metropolis()]
#' @references Diaconis, P. and B. Sturmfels (1998). Algebraic Algorithms for
#' Sampling from Conditional Distributions. \emph{The Annals of Statistics}
#' 26(1), pp.363-397.
#' @references Drton, M., B. Sturmfels, and S. Sullivant (2009). \emph{Lectures
#' on Algebraic Statistics}, Basel: Birkhauser Verlag AG.
#' @references Aoki, S., H. Hara, and A. Takemura (2012). \emph{Markov Bases in
#' Algebraic Statistics}, Springer.
#' @references Agresti, A. (2002). \emph{Categorical Data Analysis}, Basel: John
#' Wiley & Sons, 2ed.
#' @references Agresti, A. (1992). A Survey of Exact Inference for Contingency
#' Tables \emph{Statistical Science} 7(1), pp.131-153.
#' @references Read, T. and Cressie, N. (1998). \emph{Goodness-of-Fit Statistics
#' for Discrete Multivariate Data}, Springer-Verlag.
#' @examples
#'
#' \dontrun{ requires LattE and 4ti2
#'
#'
#' ## handedness introductory example
#' ############################################################
#'
#' data(handy); handy
#'
#' (out <- loglinear(~ Gender + Handedness, data = handy))
#'
#' # you can also specify the same model using variable indices...
#' (out <- loglinear(~ 1 + 2, data = handy))
#'
#' # ... or as a list of facets given by indices
#' (out <- loglinear(list(1, 2), data = handy))
#'
#' # ... or as a list of facets given by name
#' (out <- loglinear(list("Gender", "Handedness"), data = handy))
#'
#' # ... and even via a pre-computed configuration matrix
#' # this method does come with somewhat reduced output
#' A <- hmat(c(2, 2), 1:2)
#' (out <- loglinear(A, data = handy))
#'
#'
#'
#' # loglinear performs the same tasks as loglin and loglm,
#' # but loglinear gives the exact test p values and more goodness-of-fit statistics
#' stats::loglin(handy, list(1, 2))
#' MASS::loglm(~ Gender + Handedness, data = handy)
#' # loglm is just a wrapper of loglin
#'
#' # we can check loglinear's output with
#' fisher.test(handy)$p.value
#' out$p.value
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' # comparisons between loglinear, stats::loglin, and MASS::loglm
#' ############################################################
#'
#' (loglinearFit <- loglinear(~ Gender + Handedness, data = handy))
#' (loglinFit <- stats::loglin(handy, list(1, 2), fit = TRUE, param = TRUE))
#' (loglmFit <- MASS::loglm(~ Gender + Handedness, data = handy))
#'
#'
#'
#' # the expected table given the sufficient statistics can be computed
#' # via two methods, iterative proportional fitting, and the mcmc itself:
#' loglinearFit$exp # ipf
#' loglinear(~ Gender + Handedness, data = handy, method = "mcmc")$exp
#' loglinFit$fit # the equivalent in loglin; this is used by default in loglinear
#'
#'
#'
#'
#' # the parameter values of the loglinear model can be accessed
#' loglinearFit$param
#' loglinFit$param
#'
#'
#'
#'
#' # the p-value for the goodness-of-fit of the overall model is available as well :
#' # loglinear gives the exact conditional p-value
#' # (conditional on the sufficient statistics)
#' # the five numbers correspond the probability of observering a table that is
#' # "more weird" than the observed table, where "more weird" is determined
#' # by having a larger X2 value (or G2, FT, CR, or NM)
#' loglinearFit$p.value
#'
#' # in this case (a 2x2 table with the independence model), we can check that
#' # the above p-values are coorect up to Monte Carlo error
#' fisher.test(handy)$p.value
#'
#'
#' # loglin gives the p-values using the unconditional asymptotic distribution:
#' # note that the two are quite different in this case, although the conclusion
#' # is the same
#' c(
#' "X2" = pchisq(loglinFit$pearson, df = loglinFit$df, lower.tail = FALSE),
#' "G2" = pchisq(loglinFit$lrt, df = loglinFit$df, lower.tail = FALSE)
#' )
#'
#' # mid p-values are available as well:
#' loglinearFit$mid.p.value # the mid (exact conditional) p-value is also available
#'
#'
#'
#'
#' # the test statistics based on the observed table and the expected
#' # table under the model are available
#' loglinearFit$statistic
#' c(X2 = loglinFit$pearson, G2 = loglinFit$lrt) # loglin only gives X2 and G2
#' # note that PR is un-normalized log-probability
#'
#'
#'
#' # the markov moves used for the proposal distribution of the metropolis-hastings
#' # algorithm are returned. the proposal distribution is uniform on +/-
#' # the moves added to the current table
#' loglinearFit$moves
#' # they are easier understood as tables
#' vec2tab(loglinearFit$moves, dim(handy))
#' # notice that the marginals stay fixed:
#' handy + vec2tab(loglinearFit$moves, dim(handy))
#'
#'
#'
#'
#' # these were computed as the markov basis of the integer matrix
#' # (for different models, different schemes may be employed)
#' loglinearFit$A
#' markov(loglinearFit$A)
#' loglinearFit$moves
#'
#'
#'
#'
#' # the moves are also sometimes written in tableau form (LAS p.13)
#' tableau(loglinearFit$moves, dim(handy))
#' # that's +1 the the table in elements [1,1] and [2,2]
#' # and -1 in the table in elements [1,2] and [2,1]
#'
#'
#'
#'
#' # the acceptance probability of the MCMC is retained
#' loglinearFit$accept_prob
#'
#'
#'
#'
#' # various model assessment measures are also available
#' loglinearFit$quality
#'
#'
#'
#'
#' # the number of independent parameters per term are in df
#' loglinearFit$df
#'
#'
#'
#'
#' # as an added help, you may find the visuals in vcd useful:
#' library(vcd)
#' mosaic(~ Gender + Handedness, data = handy, shade = TRUE, legend = TRUE)
#' # note mosaic's use of the asymptotic X^2 test
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' ## politics example - with computing the exact p value by hand
#' ############################################################
#'
#' data(politics); politics
#'
#' (out <- loglinear(~ Personality + Party, data = politics))
#' loglinFit <- stats::loglin(politics, as.list(1:2), fit = TRUE, param = TRUE)
#'
#' out$p.value
#' # exact without monte-carlo error
#' sum(dhyper(c(0:3,6:9), 10, 10, 9))
#' fisher.test(politics)$p.value
#' round(dhyper(0:9, 10, 10, 9), 4)
#'
#' # we can sample from the hypergeometric distribution on the fiber using
#' # rhyper
#' rhyper(100, 10, 10, 9)
#'
#'
#' # comparisons :
#' out$exp
#' loglinFit$fit
#'
#' out$param
#' loglinFit$param
#'
#' out$p.value # exact
#' c(
#' X2 = pchisq(loglinFit$pearson, df = loglinFit$df, lower.tail = FALSE),
#' G2 = pchisq(loglinFit$lrt, df = loglinFit$df, lower.tail = FALSE)
#' ) # asymptotic approximation
#' fisher.test(politics)$p.value # the exact conditional p-value
#'
#' out$statistic # accurate to monte carlo error
#' c(X2 = loglinFit$pearson, G2 = loglinFit$lrt)
#'
#' vcd::mosaic(~ Personality + Party, data = politics, shade = TRUE, legend = TRUE)
#'
#'
#'
#' # alternative model specifications :
#' loglinear(~ Personality + Party, data = politics)
#' loglinear(~ 1 + 2, data = politics)
#' loglinear(list(1, 2), data = politics)
#' loglinear(list("Personality", "Party"), data = politics)
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' ## eyeHairColor from the Diaconis and Sturmfels reference
#' ############################################################
#'
#' data(HairEyeColor)
#' eyeHairColor <- margin.table(HairEyeColor, 2:1)
#'
#' out <- loglinear(~ Eye + Hair, data = eyeHairColor)
#'
#' # the default fisher.test doesn't work even with workspace = 2E9
#' # (with over 4.5Gb in memory) because it is trying to enumerate the fiber.
#' #fisher.test(eyeHairColor, workspace = 2E9)
#' # it can, however, compute Monte Carlo p-values for RxC tables, like loglinear
#' fisher.test(eyeHairColor, simulate.p.value = TRUE, B = 1e6)
#' out$p.value
#'
#' # you can see the markov moves used in the mcmc in tableau notation
#' tableau(out$moves, dim(eyeHairColor))
#'
#' # library(vcd)
#' # mosaic(~ Eye + Hair, data = HairEyeColor, shade = TRUE, legend = TRUE)
#'
#'
#'
#'
#'
#'
#' ## abortion preference example from the
#' ## Diaconis and Sturmfels reference pp. 379--381
#' ## a no 3-way interaction model
#' ############################################################
#'
#' data(abortion); abortion
#'
#' (loglinearFit <- loglinear(subsets(1:3, 2), data = abortion,
#' iter = 10000, burn = 50000, thin = 50
#' ))
#' loglinFit <- loglin(abortion, subsets(1:3, 2), fit = TRUE, param = TRUE)
#'
#' vec2tab(rowMeans(loglinearFit$steps), dim(abortion)) # cf. p. 380
#' loglinFit$fit
#'
#' all.equal(loglinearFit$param, loglinFit$param)
#'
#'
#'
#' qqplot(rchisq(1055, df = 8), out$sampsStats$X2s)
#' curve(1*x, from = 0, to = 30, add = TRUE, col = "red")
#'
#' ( nMoves <- 2*ncol(out$moves) ) # DS uses 110
#' # (the markov basis is larger than it needs to be)
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' ## loglin no three-way interaction model example
#' ############################################################
#'
#' # the help for fits the no three-way interaction model on HairEyeColor,
#' # finds a .66196 p-value using the asymptotic distribution, and concludes
#' # a good fit:
#' data(HairEyeColor)
#'
#' loglinearFit <- loglinear(subsets(1:3, 2), data = HairEyeColor)
#' loglinFit <- loglin(HairEyeColor, subsets(1:3, 2), fit = TRUE, param = TRUE)
#'
#'
#'
#'
#' # p values
#' loglinearFit$p.value
#' pchisq(loglinFit$lrt, loglinFit$df, lower.tail = FALSE) # see ?loglin
#'
#' # test statistics
#' loglinearFit$statistic
#' c(X2 = loglinFit$pearson, G2 = loglinFit$lrt)
#'
#' # fits (estimated tables)
#' loglinearFit$obs
#' round(loglinearFit$exp, 1)
#' round(loglinFit$fit, 1)
#'
#'
#' # checking the autocorrelation of mcmc
#' acf(loglinearFit$sampsStats$PRs)
#'
#' # poor mixing is a known limitation of markov bases strategies
#' # one strategy is to try to thin the mcmc
#' loglinearFit <- loglinear(subsets(1:3, 2), data = HairEyeColor, thin = 100)
#' acf(loglinearFit$sampsStats$PRs) # got it! (overkill, actually)
#'
#'
#' # the slight differences in loglinFit$fit and loglinearFit$exp (both done with ipf from loglin)
#' # are due to differences in variable order:
#' loglin(HairEyeColor, subsets(1:3, 2), fit = TRUE)$fit
#' loglin(HairEyeColor, subsets(1:3, 2)[c(1,3,2)], fit = TRUE)$fit
#'
#' # let's look at a few model moves
#' vec2tab(loglinearFit$moves[,1], dim(HairEyeColor))
#' vec2tab(loglinearFit$moves[,50], dim(HairEyeColor))
#' -vec2tab(loglinearFit$moves[,50], dim(HairEyeColor))
#'
#' # they contribute 0 to the marginals of the table
#' # (the sufficient statistics of the model)
#' exampleMove <- loglinearFit$moves[,50]
#' vec2tab(exampleMove, dim(HairEyeColor))
#' loglinearFit$A %*% exampleMove
#'
#' # two tables with same sufficient statistics
#' HairEyeColor
#' HairEyeColor + vec2tab(exampleMove, dim(HairEyeColor))
#'
#' # here are the sufficient statistics:
#' loglinearFit$A %*% tab2vec(HairEyeColor)
#' loglinearFit$A %*% tab2vec(HairEyeColor + vec2tab(exampleMove, dim(HairEyeColor)))
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' ## a table with positive marginals but no MLE for
#' ## the no-three way interaction model
#' ############################################################
#'
#'
#' data(haberman); haberman
#'
#' mod <- loglinear(subsets(1:3, 2), data = haberman)
#'
#' loglinFit <- loglin(haberman, subsets(1:3, 2), param = TRUE, fit = TRUE)
#' loglinFit$fit
#' loglinFit$param
#' c(X2 = loglinFit$pearson, G2 = loglinFit$lrt)
#'
#' loglinearFit <- loglinear(subsets(1:3, 2), data = haberman, method = "mcmc")
#' loglinearFit$exp
#' loglinearFit$param
#' loglinearFit$statistic
#'
#' A <- hmat(rep(2, 3), subsets(1:3, 2))
#' count_tables(haberman, A) # there's only one table in the fiber!
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' ## an example from agresti, p.322
#' ############################################################
#'
#' data(drugs); drugs
#' ftable(aperm(drugs, c(3, 1, 2))) # = table 8.3
#'
#' out <- loglinear(~ Alcohol + Cigarette + Marijuana, data = drugs)
#' matrix(round(aperm(out$exp, c(2,1,3)), 1), byrow = FALSE)
#'
#' loglin(drugs, as.list(1:3), fit = TRUE)$fit
#' loglin(drugs, as.list(1:3), param = TRUE)$param
#'
#' # # the saturated model issues a warning from markov, but works :
#' # out <- loglinear(~ Alcohol * Cigarette * Marijuana, data = drugs)
#' # matrix(round(aperm(out$exp, c(2,1,3)), 1), byrow = FALSE) # = the data
#'
#'
#' ftable(aperm(out$exp, c(3,1,2)))
#'
#' stats <- loglin(drugs, as.list(1:3), fit = TRUE, param = TRUE)
#'
#'
#' ## considered via glm
#'
#' df <- as.data.frame(drugs)
#' mod <- glm(Freq ~ Alcohol + Cigarette + Marijuana, data = df, family = poisson)
#' summary(mod)
#' mod$fitted.values
#'
#'
#' # the same can be done with glm :
#'
#' mod <- glm(
#' Freq ~ Alcohol + Cigarette + Marijuana,
#' data = as.data.frame(drugs), family = poisson
#' )
#' summary(mod)
#' matrix(round(mod$fitted.values[c(1,3,2,4,5,7,6,8)],1))
#'
#'
#'
#' mod <- glm(
#' Freq ~ Alcohol * Cigarette + Marijuana,
#' data = as.data.frame(drugs), family = poisson
#' )
#' summary(mod)
#' matrix(round(mod$fitted.values[c(1,3,2,4,5,7,6,8)],1))
#'
#'
#' mod <- glm(
#' Freq ~ Alcohol * Cigarette * Marijuana,
#' data = as.data.frame(drugs), family = poisson
#' )
#' summary(mod)
#' matrix(round(mod$fitted.values[c(1,3,2,4,5,7,6,8)],1))
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' }
#'
#'
loglinear <- function(model, data,
init = tab2vec(data),
iter = 1e4, burn = 1000,
thin = 10, engine = c("C++","R"),
method = c("ipf", "mcmc"), moves,
...)
{
## set/check args
##################################################
engine <- match.arg(engine)
method <- match.arg(method)
argList <- as.list(match.call(expand.dots = TRUE))[-1]
if("formula" %in% names(argList)){
.Deprecated(msg =
'the formula argument is deprecated, please use "model" instead.'
)
}
## reshape data
##################################################
data <- suppressMessages(teshape(data, "tab"))
p <- length(dim(data))
nCells <- length(data)
## if a pure array is given, give names for later
if(is.array(data) && is.null(dimnames(data))) data <- array2tab(data)
## other basic objects
varsNlevels <- dimnames(data)
vars <- names(varsNlevels)
## check for sampling zeros
##################################################
if(any(data == 0L)) message(
"Care ought be taken with tables with sampling zeros to ensure the MLE exists."
)
## parse model specification (formula for vector of r_k's)
##################################################
modelGivenByMatrix <- ifelse(is.matrix(model), TRUE, FALSE)
if(modelGivenByMatrix){
A <- model
} else {
# if it's a formula ~ Gender + Handedness or ~ 1 + 2, convert to list
if(is.formula(model)){
## parse formula
fString <- as.character(model)
predString <- fString[2]
## make list of facets
model <- strsplit(predString, " \\+ ")[[1]]
model <- strsplit(model, " \\* ")
}
# make facets (list of index vecs); if model specified with variable
# names, convert them to indices
if(all(unlist(model) %in% vars)){ # variable names
varname2index <- 1:p
names(varname2index) <- vars
facets <- lapply(model, function(varsInFacet) varname2index[varsInFacet])
} else if(all(unlist(model) %in% 1:length(vars))){ # by indices
facets <- lapply(model, as.integer) # to fix the ~ 1 + 2 case, parsed as chars
} else {
stop("Invalid model specification, see ?loglinear")
}
# make configuration (model) matrix
A <- hmat(dim(data), facets)
}
## construct A matrix and compute moves
##################################################
if(missing(moves) && has_4ti2()){
message("Computing Markov moves (4ti2)... ", appendLF = FALSE)
moves <- markov(A, p = "arb")
message("done.", appendLF = TRUE)
} else if(missing(moves) && has_4ti2()){
warning(
"No moves were provided and has_4ti2() = FALSE.\n",
" SIS moves will be used; estimates will likely be biased.\n",
" Consider using rmove() to generate SIS moves in advance.",
immediate. = TRUE
)
message("Computing 1000 SIS moves... ", appendLF = FALSE)
moves <- rmove(n = 1000, A = A, b = A %*% tab2vec(data), ...)
message("done.", appendLF = TRUE)
} else if(is.character(moves)){
movesMat <- NULL
stopifnot(all(moves %in% c("lattice", "markov", "groebner", "grobner", "graver", "sis")))
if("lattice" %in% moves) movesMat <- cbind(movesMat, zbasis(A, p = "arb"))
if("markov" %in% moves) movesMat <- cbind(movesMat, markov(A, p = "arb"))
if("groebner" %in% moves) movesMat <- cbind(movesMat, groebner(A, p = "arb"))
if("grobner" %in% moves) movesMat <- cbind(movesMat, groebner(A, p = "arb"))
if("graver" %in% moves) stop("graver not yet implemented.")
moves <- movesMat
}
stopifnot(is.array(moves))
## run metropolis-hastings
##################################################
init <- unname(init) # init
out <- metropolis(init, moves, iter = iter, burn = burn, thin = thin, engine = engine)
## compute data chi square
##################################################
if(modelGivenByMatrix && method == "ipf"){
message(
"Iterative proportional fitting is not yet implemented\n",
" for models specified by configuration matrices.\n",
" Changing to method = \"mcmc\"..."
)
method <- "mcmc"
}
if(method == "ipf"){
exp <- loglin(data, facets, fit = TRUE, print = FALSE)$fit
} else if(method == "mcmc"){
exp <- vec2tab(rowMeans(out$steps), dim(data))
dimnames(exp) <- dimnames(data)
}
e <- unname(tab2vec(exp))
u <- t(t(unname(tab2vec(data))))
PR <- computeUProbsCpp(matrix(u)) # unnormd prob; numers LAS 1.1.10
X2 <- computeX2sCpp(u, e)
G2 <- computeG2sCpp(u, e)
FT <- computeCRsCpp(u, e, -.5)
CR <- computeCRsCpp(u, e, 2/3)
NM <- computeNMsCpp(u, e)
## compute MCMC chi squares
##################################################
PRs <- computeUProbsCpp(out$steps) # unnormd probs; numers LAS 1.1.10
X2s <- computeX2sCpp(out$steps, e)
G2s <- computeG2sCpp(out$steps, e)
FTs <- computeCRsCpp(out$steps, e, -.5)
CRs <- computeCRsCpp(out$steps, e, 2/3)
NMs <- computeNMsCpp(out$steps, e)
## compute parameters
##################################################
if(!modelGivenByMatrix){
# in principle, there should be one parameter for every cell.
# there are prod(dim(data)) cells.
# a good reference is BFH, p. 35 (and to a lesser extent 43)
# the prod(dim(data)[terms[[j]]] - 1) line below is like
# (I - 1) (J - 1) (K - 1)
# CDA p.79 also helpful
dimSatModel <- nCells - 1
degFreedom <- rep.int(0, 2^p) # there are 2^p possible subsets of vars, and
# therefore there are 2^p possible terms
# possibleTerms are more "types of terms" as opposed to individual terms
# for example, an entry c(1,3) would refer to all combinations of levels
# of variables 1 and 3; ie (# var 1 levels - 1) * (# var 3 levels - 1)
# individual terms (parameters)
possibleTerms <- subsets(p, include_null = TRUE)
names(possibleTerms) <- sapply(possibleTerms, paste, collapse = " ")
names(possibleTerms)[which(names(possibleTerms) == "")] <- "(Intercept)"
nVarLvls <- dim(data)
# paramsPerTerm <- lapply(possibleTerms, function(x){
# if(length(x) == 0) return(1L)
# prod(nVarLvls[x] - 1)
# })
# similarly, there are the terms in the model
termsInModel <- unique(unlist(lapply(
lapply(facets, as.character), # to avoid subsets(2)
subsets, include_null = TRUE),
recursive = FALSE
))
termsInModel <- lapply(termsInModel, as.integer)
names(termsInModel) <- sapply(termsInModel, paste, collapse = " ")
names(termsInModel)[which(names(termsInModel) == "")] <- "(Intercept)"
paramsPerTermInModel <- lapply(termsInModel, function(x){
if(length(x) == 0) return(1L)
prod(nVarLvls[x] - 1)
})
names(paramsPerTermInModel) <- unname(sapply(termsInModel, function(x){
if(length(x) == 0) return("(Intercept)")
paste(names(dimnames(data))[x], collapse = ".")
}))
nParamsInModel <- sum(unlist(paramsPerTermInModel))
dimModel <- nParamsInModel - 1 # the - 1 accounts for the overall mean
overallAsymptoticDegFreedom <- (dimSatModel - dimModel)
# compute the parameters
log_fit <- exp
log_fit[exp > 0] <- log(exp[exp > 0])
param <- as.list(rep(NA, length(termsInModel)))
names(param) <- names(paramsPerTermInModel)
for(k in seq_along(param)){
if(length(termsInModel[[k]]) == 0){
param[[k]] <- mean(log_fit)
log_fit <- log_fit - param[[k]]
} else {
param[[k]] <- apply(log_fit, termsInModel[[k]], mean)
log_fit <- sweep(log_fit, termsInModel[[k]], param[[k]])
}
}
# for every step, fit mle
# then decompose mle
# problem : they all have the same marginals, so the same
# mles!
# idea 1 : sample from the multinomial with the same sample
# size (so different marginals), estimate, then decompose
# idea 2 : bootstrap sample from the table, estimate, decompose
# i think i like idea 2 better.
# reorder the param estimates in the order of subsets
# so you have the intercept, then all first order terms, and so on
goodOrder <- sapply(
c("(Intercept)", subsets(names(dimnames(data)))),
paste, collapse = "."
)
param <- param[goodOrder[goodOrder %in% names(param)]]
out$param <- param
}
## compute residuals and model selection, agresti p.81, 216, 324
##################################################
out$residuals <- exp
out$residuals[exp > 0] <-
(data[exp > 0] - exp[exp > 0]) / sqrt(exp[exp > 0])
if(!modelGivenByMatrix){
k <- nParamsInModel # = number of params
n <- sum(data) # = sample size
L <- dmultinom(u, sum(u), e, TRUE) # maximized log-likelihood
BIC <- log(n)*k - 2*L
AIC <- 2*k - 2*L
AICc <- AIC + 2*k*(k+1)/(n-k-1)
out$df <- paramsPerTermInModel
out$quality <- c(AIC = AIC, AICc = AICc, BIC = BIC)
}
## add A matrix, p.value and return
##################################################
out$call <- match.call()
out$obs <- data
out$exp <- exp
out$A <- A
out$p.value <- c(
PR = mean(PRs <= PR),
X2 = mean(X2s >= X2),
G2 = mean(G2s >= G2),
FT = mean(FTs >= FT),
CR = mean(CRs >= CR),
NM = mean(NMs >= NM)
)
out$p.value.std.err <- c(
PR = sqrt(mean(PRs <= PR)*(1-mean(PRs <= PR))/iter),
X2 = sqrt(mean(X2s >= X2)*(1-mean(X2s >= X2))/iter),
G2 = sqrt(mean(G2s >= G2)*(1-mean(G2s >= G2))/iter),
FT = sqrt(mean(FTs >= FT)*(1-mean(FTs >= FT))/iter),
CR = sqrt(mean(CRs >= CR)*(1-mean(CRs >= CR))/iter),
NM = sqrt(mean(NMs >= NM)*(1-mean(NMs >= NM))/iter)
)
out$mid.p.value <- c(
PR = mean(PRs < PR) + mean(PRs == PR)/2,
X2 = mean(X2s > X2) + mean(X2s == X2)/2,
G2 = mean(G2s > G2) + mean(G2s == G2)/2,
FT = mean(FTs > FT) + mean(FTs == FT)/2,
CR = mean(CRs > CR) + mean(CRs == CR)/2,
NM = mean(NMs > NM) + mean(NMs == NM)/2
)
out$iter <- iter
out$burn <- burn
out$thin <- thin
out$statistic <- c(PR = PR, X2 = X2, G2 = G2, FT = FT, CR = CR, NM = NM)
out$sampsStats <- list(PRs = PRs, X2s = X2s, G2s = G2s, FTs = FTs, CRs = CRs, NMs = NMs)
out$cells <- nCells
out$method <- method
class(out) <- "loglinear"
out
}
#' @export
hierarchical <- function(...){
.Deprecated(msg = "hierarchical is deprecated, use loglinear instead.")
loglinear(...)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.