Add observation-level sample weights to poLCA

R/poLCA.R

@@ -1,7 +1,27 @@
 poLCA <-
 function(formula,data,nclass=2,maxiter=1000,graphs=FALSE,tol=1e-10,
-                na.rm=TRUE,probs.start=NULL,nrep=1,verbose=TRUE,calc.se=TRUE) {
+                na.rm=TRUE,probs.start=NULL,nrep=1,verbose=TRUE,calc.se=TRUE,
+                weights=NULL) {
+    # weights: optional observation-level sample weights (e.g. survey weights).
+    #   Either a character string naming a column of `data`, or a numeric
+    #   vector of positive weights with length nrow(data). When supplied the
+    #   model is estimated by weighted pseudo-maximum-likelihood: the E-step
+    #   posteriors are unchanged, while the M-step response probabilities,
+    #   mixing proportions / covariate coefficients and the log-likelihood
+    #   weight each observation's contribution by its weight. Standard errors
+    #   use a pseudo-ML sandwich estimator. weights=NULL (or all-1 weights)
+    #   reproduces the original unweighted estimator exactly.
     starttime <- Sys.time()
+    # resolve weights before any row subsetting
+    if (!is.null(weights)) {
+        if (is.character(weights) & (length(weights)==1)) {
+            if (!(weights %in% colnames(data))) stop("weights column not found in data")
+            weights <- data[[weights]]
+        }
+        if (length(weights) != nrow(data)) stop("weights must have length nrow(data)")
+        if (any(!is.finite(weights)) | any(weights<=0)) stop("weights must be finite and strictly positive")
+        weights <- as.numeric(weights)
+    }
     mframe <- model.frame(formula,data,na.action=NULL)
     mf <- model.response(mframe)
     if (any(mf<1,na.rm=TRUE) | any(round(mf) != mf,na.rm=TRUE)) {
@@ -11,10 +31,16 @@ function(formula,data,nclass=2,maxiter=1000,graphs=FALSE,tol=1e-10,
     outcome categories for each variable. \n\n")
         ret <- NULL
     } else {
-    data <- data[rowSums(is.na(model.matrix(formula,mframe)))==0,]
+    keeprows <- rowSums(is.na(model.matrix(formula,mframe)))==0
+    data <- data[keeprows,]
+    if (!is.null(weights)) weights <- weights[keeprows]
     if (na.rm) {
         mframe <- model.frame(formula,data)
         y <- model.response(mframe)
+        if (!is.null(weights)) {
+            naact <- attr(mframe,"na.action")
+            if (!is.null(naact)) weights <- weights[-naact]
+        }
     } else {
         mframe <- model.frame(formula,data,na.action=NULL)
         y <- model.response(mframe)
@@ -32,21 +58,25 @@ function(formula,data,nclass=2,maxiter=1000,graphs=FALSE,tol=1e-10,
     R <- nclass
     S <- ncol(x)
     if (S>1) { calc.se <- TRUE }
+    # unified internal weight vector; usew flags the weighted estimator
+    usew <- !is.null(weights) && !all(weights==1)
+    w <- if (is.null(weights)) rep(1,N) else weights
+    sumw <- sum(w)
     eflag <- FALSE
     probs.start.ok <- TRUE
     ret <- list()
     if (R==1) {
         ret$probs <- list()
         for (j in 1:J) {
             ret$probs[[j]] <- matrix(NA,nrow=1,ncol=K.j[j])
-            for (k in 1:K.j[j]) { ret$probs[[j]][k] <- sum(y[,j]==k)/sum(y[,j]>0) }
+            for (k in 1:K.j[j]) { ret$probs[[j]][k] <- sum(w*(y[,j]==k))/sum(w*(y[,j]>0)) }
         }
         ret$probs.start <- ret$probs
         ret$P <- 1
         ret$posterior <- ret$predclass <- prior <- matrix(1,nrow=N,ncol=1)
-        ret$llik <- sum(log(poLCA.ylik.C(poLCA.vectorize(ret$probs),y)) - log(.Machine$double.xmax))
+        ret$llik <- sum(w*(log(poLCA.ylik.C(poLCA.vectorize(ret$probs),y)) - log(.Machine$double.xmax)))
         if (calc.se) {
-            se <- poLCA.se(y,x,ret$probs,prior,ret$posterior)
+            se <- poLCA.se(y,x,ret$probs,prior,ret$posterior,w=w)
             ret$probs.se <- se$probs           # standard errors of class-conditional response probabilities
             ret$P.se <- se$P                   # standard errors of class population shares
         } else {
@@ -84,9 +114,9 @@ function(formula,data,nclass=2,maxiter=1000,graphs=FALSE,tol=1e-10,
                 prior <- poLCA.updatePrior(b,x,R)
                 if ((!probs.start.ok) | (is.null(probs.start)) | (!firstrun) | (repl>1)) { # only use the specified probs.start in the first nrep
                     probs <- list()
-                    for (j in 1:J) { 
+                    for (j in 1:J) {
                         probs[[j]] <- matrix(runif(R*K.j[j]),nrow=R,ncol=K.j[j])
-                        probs[[j]] <- probs[[j]]/rowSums(probs[[j]]) 
+                        probs[[j]] <- probs[[j]]/rowSums(probs[[j]])
                     }
                     probs.init <- probs
                 }
@@ -98,25 +128,36 @@ function(formula,data,nclass=2,maxiter=1000,graphs=FALSE,tol=1e-10,
                 while ((iter <= maxiter) & (dll > tol) & (!error)) {
                     iter <- iter+1
                     rgivy <- poLCA.postClass.C(prior,vp,y)      # calculate posterior
-                    vp$vecprobs <- poLCA.probHat.C(rgivy,y,vp)  # update probs
+                    if (usew) {
+                        # weighted M-step: probhat normalizes by class totals of the
+                        # posterior, so row-scaling the posterior by w yields the
+                        # weighted estimator of the response probabilities
+                        vp$vecprobs <- poLCA.probHat.C(w*rgivy,y,vp)
+                    } else {
+                        vp$vecprobs <- poLCA.probHat.C(rgivy,y,vp)  # update probs
+                    }
                     if (S>1) {
-                        dd <- poLCA.dLL2dBeta.C(rgivy,prior,x)
+                        if (usew) {
+                            dd <- poLCA.dLL2dBeta.w(rgivy,prior,x,w)
+                        } else {
+                            dd <- poLCA.dLL2dBeta.C(rgivy,prior,x)
+                        }
                         b <- b + ginv(-dd$hess) %*% dd$grad     # update betas
                         prior <- poLCA.updatePrior(b,x,R)       # update prior
                     } else {
-                        prior <- matrix(colMeans(rgivy),nrow=N,ncol=R,byrow=TRUE)
+                        prior <- matrix(colSums(w*rgivy)/sumw,nrow=N,ncol=R,byrow=TRUE)
                     }
-                    llik[iter] <- sum(log(rowSums(prior*poLCA.ylik.C(vp,y))) - log(.Machine$double.xmax))
+                    llik[iter] <- sum(w*(log(rowSums(prior*poLCA.ylik.C(vp,y))) - log(.Machine$double.xmax)))
                     dll <- llik[iter]-llik[iter-1]
                     if (is.na(dll)) {
                         error <- TRUE
                     } else if ((S>1) & (dll < -1e-7)) {
                         error <- TRUE
                     }
                 }
-                if (!error) { 
+                if (!error) {
                     if (calc.se) {
-                        se <- poLCA.se(y,x,poLCA.unvectorize(vp),prior,rgivy)
+                        se <- poLCA.se(y,x,poLCA.unvectorize(vp),prior,rgivy,w=w)
                     } else {
                         se <- list(probs=NA,P=NA,b=NA,var.b=NA)
                     }
@@ -134,7 +175,7 @@ function(formula,data,nclass=2,maxiter=1000,graphs=FALSE,tol=1e-10,
                 ret$P.se <- se$P                   # standard errors of class population shares
                 ret$posterior <- rgivy             # NxR matrix of posterior class membership probabilities
                 ret$predclass <- apply(ret$posterior,1,which.max)   # Nx1 vector of predicted class memberships, by modal assignment
-                ret$P <- colMeans(ret$posterior)   # estimated class population shares
+                ret$P <- colSums(w*ret$posterior)/sumw   # estimated class population shares (weighted when weights supplied)
                 ret$numiter <- iter-1              # number of iterations until reaching convergence
                 ret$probs.start.ok <- probs.start.ok # if starting probs specified, logical indicating proper entry format
                 if (S>1) {
@@ -159,24 +200,26 @@ function(formula,data,nclass=2,maxiter=1000,graphs=FALSE,tol=1e-10,
     ret$npar <- (R*sum(K.j-1)) + (R-1)                  # number of degrees of freedom used by the model (number of estimated parameters)
     if (S>1) { ret$npar <- ret$npar + (S*(R-1)) - (R-1) }
     ret$aic <- (-2 * ret$llik) + (2 * ret$npar)         # Akaike Information Criterion
-    ret$bic <- (-2 * ret$llik) + (log(N) * ret$npar)    # Schwarz-Bayesian Information Criterion
+    ret$bic <- (-2 * ret$llik) + (log(N) * ret$npar)    # Schwarz-Bayesian Information Criterion (uses N observations, not sum of weights)
     ret$Nobs <- sum(rowSums(y==0)==0)                   # number of fully observed cases (if na.rm=F)
     if (all(rowSums(y==0)>0)) { # if no rows are fully observed
         ret$Chisq <- NA
         ret$Gsq <- NA
         ret$predcell <- NA
     } else {
-        compy <- poLCA.compress(y[(rowSums(y==0)==0),])
+        fullobs <- rowSums(y==0)==0
+        compy <- poLCA.compress(y[fullobs,],w=w[fullobs])
         datacell <- compy$datamat
         rownames(datacell) <- NULL
-        freq <- compy$freq
-	ylik <- poLCA.ylik.C(poLCA.vectorize(ret$probs),datacell) 
+        freq <- compy$freq                     # weighted cell frequencies when weights supplied
+        ylik <- poLCA.ylik.C(poLCA.vectorize(ret$probs),datacell)
         if (!na.rm) {
-            fit <- matrix(ret$Nobs/.Machine$double.xmax * (ylik  %*% ret$P))
-            ret$Chisq <- sum((freq-fit)^2/fit) + (ret$Nobs-sum(fit)) # Pearson Chi-square goodness of fit statistic for fitted vs. observed multiway tables
+            sumw.obs <- sum(w[fullobs])
+            fit <- matrix(sumw.obs/.Machine$double.xmax * (ylik  %*% ret$P))
+            ret$Chisq <- sum((freq-fit)^2/fit) + (sumw.obs-sum(fit)) # Pearson Chi-square goodness of fit statistic for fitted vs. observed multiway tables
         } else {
-            fit <- matrix(N/.Machine$double.xmax * (ylik %*% ret$P))
-            ret$Chisq <- sum((freq-fit)^2/fit) + (N-sum(fit))
+            fit <- matrix(sumw/.Machine$double.xmax * (ylik %*% ret$P))
+            ret$Chisq <- sum((freq-fit)^2/fit) + (sumw-sum(fit))
         }
         ret$predcell <- data.frame(datacell,observed=freq,expected=round(fit,3)) # Table that gives observed vs. predicted cell counts
         ret$Gsq <- 2 * sum(freq*log(freq/fit))  # Likelihood ratio/deviance statistic
@@ -195,6 +238,7 @@ function(formula,data,nclass=2,maxiter=1000,graphs=FALSE,tol=1e-10,
         }
     }
     ret$N <- N                         # number of observations
+    ret$weights <- if (usew) w else NULL # observation-level sample weights actually used (NULL if unweighted)
     ret$maxiter <- maxiter             # maximum number of iterations specified by user
     ret$resid.df <- min(ret$N,(prod(K.j)-1))-ret$npar # number of residual degrees of freedom
     class(ret) <- "poLCA"

R/poLCA.compress.R

@@ -1,16 +1,20 @@
 poLCA.compress <-
-function(y) {
-    ym.sorted <- y[do.call(order,data.frame(y)),]
+function(y,w=NULL) {
+# w: optional observation weights; freq becomes the sum of weights in each
+# cell (equals the original cell counts when w is NULL or all 1).
+    ord <- do.call(order,data.frame(y))
+    ym.sorted <- y[ord,]
+    w.sorted <- if (is.null(w)) rep(1,nrow(ym.sorted)) else w[ord]
     vars <- ncol(ym.sorted)
     datamat <- ym.sorted[1,]
-    freq <- 1
+    freq <- w.sorted[1]
     curpos <- 1
     for (i in 2:nrow(ym.sorted)) {
         if (sum(ym.sorted[i,] == ym.sorted[i-1,])==vars) {
-            freq[curpos] <- freq[curpos]+1
+            freq[curpos] <- freq[curpos]+w.sorted[i]
         } else {
             datamat <- rbind(datamat,ym.sorted[i,])
-            freq <- c(freq,1)
+            freq <- c(freq,w.sorted[i])
             curpos <- curpos+1
         }
     }

R/poLCA.dLL2dBeta.w.R

@@ -0,0 +1,29 @@
+poLCA.dLL2dBeta.w <-
+function(rgivy,prior,x,w) {
+# Weighted (pure R) analogue of the C function d2lldbeta2 / poLCA.dLL2dBeta.C.
+# Each observation's gradient and hessian contribution is multiplied by its
+# sample weight w[i]. With w = rep(1,N) this reproduces poLCA.dLL2dBeta.C
+# exactly (up to numerical error).
+    R <- ncol(prior)
+    S <- ncol(x)
+    crank <- S*(R-1)
+    d <- rgivy - prior
+    grad <- as.vector(vapply(2:R, function(j) colSums(w * x * d[,j]), numeric(S)))
+    hess <- matrix(0,nrow=crank,ncol=crank)
+    for (j in 2:R) {
+        for (n in 2:j) {
+            coefvec <- if (n==j) {
+                w * ( -rgivy[,j]*(1-rgivy[,j]) + prior[,j]*(1-prior[,j]) )
+            } else {
+                w * ( rgivy[,j]*rgivy[,n] - prior[,j]*prior[,n] )
+            }
+            blk <- crossprod(x, x * coefvec)
+            hess[((j-2)*S+1):((j-1)*S),((n-2)*S+1):((n-1)*S)] <- blk
+            if (n!=j) {
+                hess[((n-2)*S+1):((n-1)*S),((j-2)*S+1):((j-1)*S)] <- t(blk)
+            }
+        }
+    }
+    # poLCA.dLL2dBeta.C returns hess = -(accumulated matrix); match that here
+    return(list(grad=grad,hess=-hess))
+}

R/poLCA.se.R

@@ -1,5 +1,9 @@
 poLCA.se <-
-function(y,x,probs,prior,rgivy) {
+function(y,x,probs,prior,rgivy,w=NULL) {
+# w: optional observation-level sample weights. When supplied, the variance
+# is the pseudo-maximum-likelihood sandwich estimator
+#   VCE = A^-1 B A^-1,  A = sum_i w_i s_i s_i',  B = sum_i w_i^2 s_i s_i'
+# which reduces to the original ginv(t(s) %*% s) when all w equal 1.
     J <- ncol(y)
     R <- ncol(prior)
     K.j <- sapply(probs,ncol)
@@ -23,8 +27,15 @@ function(y,x,probs,prior,rgivy) {
     if (R>1) for (r in 2:R) { s <- cbind(s,x*ppdiff[,r]) }
 
     s[is.na(s)] <- 0      # To handle missing values
-    info <- t(s) %*% s    # Information matrix
-    VCE <- ginv(info)     # VCE matrix of log-odds response probs and covariate coefficients
+    if (is.null(w) || all(w==1)) {
+        info <- t(s) %*% s    # Information matrix
+        VCE <- ginv(info)     # VCE matrix of log-odds response probs and covariate coefficients
+    } else {
+        A <- t(s) %*% (w*s)          # weighted information matrix
+        B <- t(w*s) %*% (w*s)        # meat of the sandwich
+        Ainv <- ginv(A)
+        VCE <- Ainv %*% B %*% t(Ainv) # pseudo-ML sandwich VCE
+    }
 
     # Variance of class conditional response probs using delta fn. transformation with 
     # Jacobian a block diagonal matrix (across r) of block diagonal matrices (across j)
@@ -67,9 +78,10 @@ function(y,x,probs,prior,rgivy) {
             diag(ptp[,,n]) <- prior[n,] * (1-prior[n,])
         }
         Jac.mix <- NULL
+        wjac <- if (is.null(w)) rep(1,N) else w
         for (r in 2:R) {
             for (l in 1:ncol(x)) {
-                Jac.mix <- cbind(Jac.mix,colMeans(t(ptp[,r,]) * x[,l]))
+                Jac.mix <- cbind(Jac.mix,colSums(wjac * t(ptp[,r,]) * x[,l])/sum(wjac))
             }
         }
         VCE.mix <- Jac.mix %*% VCE.beta %*% t(Jac.mix)