diff --git a/R/poLCA.R b/R/poLCA.R index 35d3949..432235d 100644 --- a/R/poLCA.R +++ b/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,6 +58,10 @@ 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() @@ -39,14 +69,14 @@ function(formula,data,nclass=2,maxiter=1000,graphs=FALSE,tol=1e-10, 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,15 +128,26 @@ 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 @@ -114,9 +155,9 @@ function(formula,data,nclass=2,maxiter=1000,graphs=FALSE,tol=1e-10, 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" diff --git a/R/poLCA.compress.R b/R/poLCA.compress.R index 7458fad..d064141 100644 --- a/R/poLCA.compress.R +++ b/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 } } diff --git a/R/poLCA.dLL2dBeta.w.R b/R/poLCA.dLL2dBeta.w.R new file mode 100644 index 0000000..82ec1f3 --- /dev/null +++ b/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)) +} diff --git a/R/poLCA.se.R b/R/poLCA.se.R index 8802f39..efbd405 100644 --- a/R/poLCA.se.R +++ b/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)