Bvr

DESCRIPTION

@@ -14,3 +14,5 @@ URL: http://dlinzer.github.com/poLCA
 LazyLoad: yes
 NeedsCompilation: yes
 Packaged: 2016-06-04 16:53:43 UTC; jeff
+Suggests: testthat
+RoxygenNote: 6.0.1

R/poLCA.R

@@ -178,7 +178,7 @@ function(formula,data,nclass=2,maxiter=1000,graphs=FALSE,tol=1e-10,
             fit <- matrix(N/.Machine$double.xmax * (ylik %*% ret$P))
             ret$Chisq <- sum((freq-fit)^2/fit) + (N-sum(fit))
         }
-        ret$predcell <- data.frame(datacell,observed=freq,expected=round(fit,3)) # Table that gives observed vs. predicted cell counts
+        ret$predcell <- data.frame(datacell,observed=freq,expected=fit) # Table that gives observed vs. predicted cell counts
         ret$Gsq <- 2 * sum(freq*log(freq/fit))  # Likelihood ratio/deviance statistic
     }
     y[y==0] <- NA

R/poLCA.bvr.R

@@ -0,0 +1,69 @@
+#' Bivariate residuals for latent class models
+#' 
+#' Calculate the "bivariate residuals" (BVRs) between pairs of variables 
+#' in a latent class model.
+#' 
+#' This function compares the model-implied (expected) counts in the crosstables
+#' of all pairs of observed dependent variables to the observed counts. For each
+#' pair, it calculates a "chi-square" statistic,
+#' 
+#' \deqn{\text{BVR} = \sum_{j, j'} \frac{(n_{jj'} - e_{jj'})^2}{e_{jj'}}},
+#' 
+#' where \eqn{n_{jj'}} are the observed counts for categories \eqn{j} and \eqn{j'} 
+#' of the variables being crosstabulated, and \eqn{e_{jj'}} are
+#' the expected counts under the latent class model. 
+#' 
+#' Note that the BVR does not follow an asymptotic chi-square distribution and
+#' for accurate p-values, parametric bootstrapping is necessary (Oberski et al. 2013).
+#' 
+#' @param fit A poLCA fit object
+#' @param tol Optional: tolerance for small expected counts
+#' @param rescale_to_df Optional: whether to divide the pairwise "chi-square" values by 
+#' the degrees of freedom of the local crosstable. Default is TRUE.
+#' @return The table of bivariate residuals
+#' @author Daniel Oberski ([email protected])
+#' @seealso \code{\link{poLCA}} for fitting the latent class model.
+#' @references 
+#' Oberski, DL, Van Kollenburg, GH and Vermunt, JK (2013). 
+#'   A Monte Carlo evaluation of three methods to detect local dependence in binary data latent class models. 
+#'   Advances in Data Analysis and Classification 7 (3), 267-279.
+#' @examples
+#' data(values)
+#' f <- cbind(A, B, C, D) ~ 1
+#' M0 <- poLCA(f,values, nclass=1, verbose = FALSE) 
+#' bvr(M0) # 12.4, 5.7, 8.3, 15.6, ... 
+bvr <- function(fit, tol = 1e-3, rescale_to_df = TRUE) {
+  stopifnot(class(fit) == "poLCA")
+
+  ov_names <- names(fit$predcell)[1:(ncol(fit$predcell) - 2)]
+  ov_combn <- combn(ov_names, 2)
+
+  get_bvr <- function(ov_pair) {
+    form_obs <- as.formula(paste0("observed ~ ", ov_pair[1], " + ", ov_pair[2]))
+    form_exp <- as.formula(paste0("expected ~ ", ov_pair[1], " + ", ov_pair[2]))
+
+    counts_obs <- xtabs(form_obs, data = fit$predcell)
+    counts_exp <- xtabs(form_exp, data = fit$predcell)
+    counts_exp <- ifelse(counts_exp < tol, tol, counts_exp) # Prevent Inf/NaN
+
+    bvr_df <- prod(dim(counts_exp) - 1)
+    bvr_value <- sum((counts_obs - counts_exp)^2 / counts_exp)
+
+    if(rescale_to_df) bvr_value <- bvr_value / bvr_df
+
+    attr(bvr_value, "df") <- bvr_df
+
+    bvr_value
+  }
+
+  bvr_pairs <- apply(ov_combn, 2, get_bvr)
+
+  attr(bvr_pairs, "rescale_to_df") <- rescale_to_df
+  attr(bvr_pairs, "class") <- "dist"
+  attr(bvr_pairs, "Size") <- length(ov_names)
+  attr(bvr_pairs, "Labels") <- ov_names
+  attr(bvr_pairs, "Diag") <- FALSE
+  attr(bvr_pairs, "Upper") <- FALSE
+
+  bvr_pairs
+}

tests/testthat.R

@@ -0,0 +1,4 @@
+library(testthat)
+library(poLCA)
+
+test_check("poLCA")

tests/testthat/test-bvr.r

@@ -0,0 +1,33 @@
+library(poLCA)
+context("BVRs for latent class models")
+
+
+test_that("BVRs for example model are as expected (no covariates)", {
+  set.seed(201708)
+
+  data(values)
+  f <- cbind(A,B,C,D)~1
+  M0 <- poLCA(f,values,nclass=1, verbose = FALSE) # log-likelihood: -543.6498
+  M1 <- poLCA(f,values,nclass=2, verbose = FALSE) # log-likelihood: -504.4677
+
+  bvr_0 <- bvr(M0)
+  bvr_1 <- bvr(M1)
+
+  expect_that(all(abs(bvr_0 - c(12.378947,  5.704402,  8.308752, 15.586100, 23.562456, 23.779559)) < 1e-6), is_true() )
+  expect_that(all(abs(bvr_1 - c(1.172161691, 0.003303647, 0.010097938, 0.037913535, 0.016503607, 0.055858509)) < 1e-6), is_true() )
+
+})
+
+test_that("BVRs for example model are as expected (with covariates)", {
+  set.seed(201708)
+
+  data(election)
+
+  f2a <- cbind(MORALG, CARESG, DISHONG, INTELG) ~ PARTY
+
+  fit_nes <-poLCA(f2a, election, nclass = 3, nrep = 5, verbose = T)  
+
+  bvr_nes <- bvr(fit_nes)
+
+  expect_that(all(abs(bvr_nes - c(5.773705,  6.799946,  7.369901,  4.684608,  5.205810, 11.022726)) < 1e-6), is_true() )
+})