|
| 1 | +### R-side limited-information (pairwise) posterior predictive check for |
| 2 | +### models with ordinal indicators. |
| 3 | +### |
| 4 | +### Motivation: the default Stan-computed ppp compares the fitted model to |
| 5 | +### a saturated model on the continuous data-augmented y* used internally |
| 6 | +### for ordinal likelihoods. Because that augmentation already respects the |
| 7 | +### fitted model's own thresholds, the comparison understates ordinal lack |
| 8 | +### of fit. This file instead compares observed vs model-implied pairwise |
| 9 | +### summaries (contingency tables for ordinal-ordinal pairs, a conditional/ |
| 10 | +### polyserial statistic for ordinal-continuous pairs, correlations for |
| 11 | +### continuous-continuous pairs), separately for the real data and for a |
| 12 | +### posterior-predictive replicate, per retained MCMC draw. |
| 13 | +### |
| 14 | +### Scope: single-level (single- or multi-group) target="stan"/"cmdstan" |
| 15 | +### models only; see R/ppp.R for the dispatcher that keeps this out of the |
| 16 | +### two-level and continuous-only paths. |
| 17 | + |
| 18 | +## classify every pair of observed variables (assumes structurally |
| 19 | +## identical ov.names/types across groups, as blavaan requires) as |
| 20 | +## oo (ordinal-ordinal), oc (ordinal-continuous), or cc (continuous-continuous) |
| 21 | +li_pairtypes <- function(lavdata, lavsamplestats) { |
| 22 | + ov.names <- lavdata@ov.names[[1]] |
| 23 | + nvar <- length(ov.names) |
| 24 | + ov.type <- lavdata@ov$type[match(ov.names, lavdata@ov$name)] |
| 25 | + is.ord <- ov.type == "ordered" |
| 26 | + |
| 27 | + ## a cc pair is dropped only if fixed.x in *every* group (conservative: |
| 28 | + ## a pair that is fixed.x in some but not all groups is kept) |
| 29 | + ngroups <- lavsamplestats@ngroups |
| 30 | + fixedx <- rep(TRUE, nvar) |
| 31 | + for (g in 1:ngroups) { |
| 32 | + x.idx <- lavsamplestats@x.idx[[g]] |
| 33 | + thisfx <- rep(FALSE, nvar) |
| 34 | + if (length(x.idx) > 0) thisfx[x.idx] <- TRUE |
| 35 | + fixedx <- fixedx & thisfx |
| 36 | + } |
| 37 | + |
| 38 | + oo <- oc <- cc <- matrix(numeric(0), 0, 2) |
| 39 | + if (nvar >= 2) { |
| 40 | + allpairs <- combn(nvar, 2) |
| 41 | + for (r in seq_len(ncol(allpairs))) { |
| 42 | + i <- allpairs[1, r]; j <- allpairs[2, r] |
| 43 | + if (is.ord[i] && is.ord[j]) { |
| 44 | + oo <- rbind(oo, c(i, j)) |
| 45 | + } else if (is.ord[i] || is.ord[j]) { |
| 46 | + ordv <- if (is.ord[i]) i else j |
| 47 | + contv <- if (is.ord[i]) j else i |
| 48 | + oc <- rbind(oc, c(ordv, contv)) |
| 49 | + } else if (!(fixedx[i] && fixedx[j])) { |
| 50 | + cc <- rbind(cc, c(i, j)) |
| 51 | + } |
| 52 | + } |
| 53 | + } |
| 54 | + |
| 55 | + list(ov.names = ov.names, is.ord = is.ord, oo = oo, oc = oc, cc = cc) |
| 56 | +} |
| 57 | + |
| 58 | +## fixed (non-draw-dependent) observed-data summaries for oc and cc pairs. |
| 59 | +## oo pairs need no bespoke summary here -- handled by lavTables() |
| 60 | +## directly off lavdata@X inside li_oo_stat(). |
| 61 | +li_obs_summary <- function(lavdata, pairs) { |
| 62 | + ngroups <- lavdata@ngroups |
| 63 | + oc_obs <- vector("list", ngroups) |
| 64 | + cc_obs <- vector("list", ngroups) |
| 65 | + for (g in 1:ngroups) { |
| 66 | + X <- lavdata@X[[g]] |
| 67 | + if (nrow(pairs$oc) > 0) { |
| 68 | + oc_obs[[g]] <- lapply(seq_len(nrow(pairs$oc)), function(r) { |
| 69 | + ordv <- pairs$oc[r, 1]; contv <- pairs$oc[r, 2] |
| 70 | + ok <- complete.cases(X[, c(ordv, contv)]) |
| 71 | + list(ord = X[ok, ordv], cont = X[ok, contv]) |
| 72 | + }) |
| 73 | + } |
| 74 | + if (nrow(pairs$cc) > 0) { |
| 75 | + cc_obs[[g]] <- lapply(seq_len(nrow(pairs$cc)), function(r) { |
| 76 | + i <- pairs$cc[r, 1]; j <- pairs$cc[r, 2] |
| 77 | + ok <- complete.cases(X[, c(i, j)]) |
| 78 | + list(r = cor(X[ok, i], X[ok, j]), n = sum(ok)) |
| 79 | + }) |
| 80 | + } |
| 81 | + } |
| 82 | + list(oc = oc_obs, cc = cc_obs) |
| 83 | +} |
| 84 | + |
| 85 | +## build a non-refit ("fake fitted") lavaan object from one posterior |
| 86 | +## draw's filled-in lavmodel and a data list (either the original observed |
| 87 | +## data, or a posterior-predictive replicate from postdata()). Mirrors the |
| 88 | +## existing hack in postpred.R's "other measures" branch (ugly hack to |
| 89 | +## avoid lav_samplestats_from_data: reconstruct data + call lavaan()). |
| 90 | +li_build_fake_lavobj <- function(datalist, lavmodel_i, lavpartable, lavoptions, |
| 91 | + lavcache, lavdata) { |
| 92 | + DATA.X <- do.call("rbind", datalist) |
| 93 | + colnames(DATA.X) <- lavdata@ov.names[[1L]] |
| 94 | + DATA <- as.data.frame(DATA.X) |
| 95 | + |
| 96 | + lavoptions2 <- lavoptions |
| 97 | + lavoptions2$verbose <- FALSE |
| 98 | + lavoptions2$estimator <- if (lavmodel_i@categorical) "DWLS" else "ML" |
| 99 | + lavoptions2$se <- "none" |
| 100 | + lavoptions2$test <- "standard" |
| 101 | + lavoptions2$optim.method <- "none" |
| 102 | + if ("control" %in% slotNames(lavmodel_i)) { |
| 103 | + lavmodel_i@control <- list(optim.method = "none") |
| 104 | + } |
| 105 | + |
| 106 | + ngroups <- length(datalist) |
| 107 | + if (ngroups > 1L) { |
| 108 | + DATA[, lavdata@group] <- rep(lavdata@group.label, |
| 109 | + times = sapply(datalist, nrow)) |
| 110 | + out <- lavaan(slotOptions = lavoptions2, slotParTable = lavpartable, |
| 111 | + slotSampleStats = NULL, slotData = NULL, |
| 112 | + slotModel = lavmodel_i, slotCache = lavcache, |
| 113 | + data = DATA, group = lavdata@group) |
| 114 | + } else { |
| 115 | + out <- lavaan(slotOptions = lavoptions2, slotParTable = lavpartable, |
| 116 | + slotSampleStats = NULL, slotData = NULL, |
| 117 | + slotModel = lavmodel_i, slotCache = lavcache, |
| 118 | + data = DATA) |
| 119 | + } |
| 120 | + ## lavaan() with optim.method="none" leaves @implied empty |
| 121 | + implied <- lav_model_implied(lavmodel_i, |
| 122 | + delta = (lavmodel_i@parameterization == "delta")) |
| 123 | + ## blavaan forces theta parameterization for ordinal models |
| 124 | + ## (R/blavaan.R), so implied$cov's diagonal is not 1 for ordinal |
| 125 | + ## variables' underlying y* -- but implied$th is already on the |
| 126 | + ## standardized (unit-variance) scale (lav_tables_pairwise_model_pi() |
| 127 | + ## uses th directly, with no rescaling of its own). Convert cov to a |
| 128 | + ## genuine correlation matrix so the two are on a consistent scale; |
| 129 | + ## otherwise lavTables()/pbivnorm sees spurious "correlations" > 1. |
| 130 | + implied$cov <- lapply(implied$cov, cov2cor) |
| 131 | + out@implied <- implied |
| 132 | + out |
| 133 | +} |
| 134 | + |
| 135 | +## thin wrapper around lavTables()'s pairwise (2-way) X2 table, |
| 136 | +## summed across all ordinal-ordinal pairs. lavTables() enumerates ordinal |
| 137 | +## pairs itself off fake_lavobj@Data@ov$type -- no pair list needed here. |
| 138 | +## A draw with an out-of-[-1,1] implied correlation (rare, but possible for |
| 139 | +## a poorly-mixed/divergent draw) makes lavTables()'s pbivnorm call error |
| 140 | +## hard rather than degrade gracefully -- treat that draw's oo contribution |
| 141 | +## as unusable (NA), mirroring get_ll_ord()'s try()-wrapped NA convention |
| 142 | +## for its own mnormt::sadmvn() calls (R/blav_model_loglik.R). |
| 143 | +li_oo_stat <- function(fake_lavobj) { |
| 144 | + tf <- try(lavTables(fake_lavobj, dimension = 2L, type = "table", |
| 145 | + statistic = "X2"), silent = TRUE) |
| 146 | + if (inherits(tf, "try-error")) return(NA_real_) |
| 147 | + if (is.null(tf) || nrow(tf) == 0) return(0) |
| 148 | + sum(tf$X2, na.rm = TRUE) |
| 149 | +} |
| 150 | + |
| 151 | +## polyserial-style discrepancy for one ordinal-continuous pair (one group). |
| 152 | +## No natural tabulation exists once one variable is continuous, so this is |
| 153 | +## a per-case deviance: -2*sum(log(conditional probability of the observed |
| 154 | +## category)), given the conditional distribution of the ordinal item's |
| 155 | +## underlying y* conditional on the paired continuous variable's value |
| 156 | +## (standard bivariate-normal conditioning). A squared-Pearson-residual |
| 157 | +## form (obs-p)^2/p was tried first but rejected: it blows up whenever any |
| 158 | +## single case's observed category has a small conditional probability |
| 159 | +## (routine in the tails, even for a correctly-specified model, once |
| 160 | +## multiplied out over many cases) -- a per-case deviance only grows |
| 161 | +## logarithmically as p -> 0, so one unusual case can't dominate the sum. |
| 162 | +li_oc_stat <- function(ord, cont, mu_o, mu_c, s_oo, s_cc, s_oc, tau_raw, |
| 163 | + floor = 1e-8) { |
| 164 | + n <- length(ord) |
| 165 | + if (n == 0) return(0) |
| 166 | + rho <- s_oc / sqrt(s_oo * s_cc) |
| 167 | + mu_cond <- mu_o + rho * sqrt(s_oo / s_cc) * (cont - mu_c) |
| 168 | + sd_cond <- sqrt(s_oo * (1 - rho^2)) |
| 169 | + cuts <- c(-Inf, tau_raw, Inf) |
| 170 | + |
| 171 | + upper <- cuts[ord + 1L] |
| 172 | + lower <- cuts[ord] |
| 173 | + catprob_obs <- pnorm((upper - mu_cond) / sd_cond) - |
| 174 | + pnorm((lower - mu_cond) / sd_cond) |
| 175 | + catprob_obs <- pmax(catprob_obs, floor) |
| 176 | + |
| 177 | + -2 * sum(log(catprob_obs)) |
| 178 | +} |
| 179 | + |
| 180 | +## continuous-continuous pair discrepancy: squared standardized residual |
| 181 | +## between the sample and model-implied correlation. |
| 182 | +li_cc_stat <- function(r_data, r_model, n) { |
| 183 | + if (n == 0 || is.na(r_data)) return(0) |
| 184 | + n * (r_data - r_model)^2 |
| 185 | +} |
| 186 | + |
| 187 | +## per-draw discrepancy, observed vs replicated data, summed across all |
| 188 | +## oo/oc/cc pairs. |
| 189 | +li_draw_stat <- function(postsamp, lavobject, lavmodel_orig, lavdata, |
| 190 | + lavsamplestats, lavpartable, lavoptions, lavcache, |
| 191 | + samp.index, chain.num, pairs, obs_summary, has_oo) { |
| 192 | + lavmodel_i <- fill_params(postsamp, lavmodel_orig, lavpartable) |
| 193 | + implied <- lav_model_implied(lavmodel_i, |
| 194 | + delta = (lavmodel_i@parameterization == "delta")) |
| 195 | + |
| 196 | + gen <- postdata(samp.indices = samp.index, chain.num = chain.num, |
| 197 | + lavmodel = lavmodel_orig, lavdata = lavdata, |
| 198 | + lavjags = lavobject@external$mcmcout, lavpartable = lavpartable, |
| 199 | + lavsamplestats = lavsamplestats, lavobject = lavobject) |
| 200 | + dataX_rep <- gen[[1]] |
| 201 | + |
| 202 | + ngroups <- lavdata@ngroups |
| 203 | + T_obs <- 0 |
| 204 | + T_rep <- 0 |
| 205 | + |
| 206 | + if (has_oo) { |
| 207 | + obs_fake <- li_build_fake_lavobj(lavdata@X, lavmodel_i, lavpartable, |
| 208 | + lavoptions, lavcache, lavdata) |
| 209 | + T_obs <- T_obs + li_oo_stat(obs_fake) |
| 210 | + |
| 211 | + rep_fake <- li_build_fake_lavobj(dataX_rep, lavmodel_i, lavpartable, |
| 212 | + lavoptions, lavcache, lavdata) |
| 213 | + T_rep <- T_rep + li_oo_stat(rep_fake) |
| 214 | + } |
| 215 | + |
| 216 | + for (g in 1:ngroups) { |
| 217 | + Xg_rep <- dataX_rep[[g]] |
| 218 | + mu <- as.numeric(implied$mean[[g]]) |
| 219 | + Sigma <- (implied$cov[[g]] + t(implied$cov[[g]])) / 2 |
| 220 | + th <- implied$th[[g]] |
| 221 | + th.idx <- lavmodel_i@th.idx[[g]] |
| 222 | + |
| 223 | + if (nrow(pairs$oc) > 0) { |
| 224 | + for (r in seq_len(nrow(pairs$oc))) { |
| 225 | + ordv <- pairs$oc[r, 1]; contv <- pairs$oc[r, 2] |
| 226 | + tau_raw <- mu[ordv] + th[th.idx == ordv] * sqrt(Sigma[ordv, ordv]) |
| 227 | + |
| 228 | + oo_obs <- obs_summary$oc[[g]][[r]] |
| 229 | + T_obs <- T_obs + li_oc_stat(oo_obs$ord, oo_obs$cont, |
| 230 | + mu[ordv], mu[contv], |
| 231 | + Sigma[ordv, ordv], Sigma[contv, contv], |
| 232 | + Sigma[ordv, contv], tau_raw) |
| 233 | + |
| 234 | + ok <- complete.cases(Xg_rep[, c(ordv, contv)]) |
| 235 | + T_rep <- T_rep + li_oc_stat(Xg_rep[ok, ordv], Xg_rep[ok, contv], |
| 236 | + mu[ordv], mu[contv], |
| 237 | + Sigma[ordv, ordv], Sigma[contv, contv], |
| 238 | + Sigma[ordv, contv], tau_raw) |
| 239 | + } |
| 240 | + } |
| 241 | + |
| 242 | + if (nrow(pairs$cc) > 0) { |
| 243 | + for (r in seq_len(nrow(pairs$cc))) { |
| 244 | + iv <- pairs$cc[r, 1]; jv <- pairs$cc[r, 2] |
| 245 | + r_model <- Sigma[iv, jv] / sqrt(Sigma[iv, iv] * Sigma[jv, jv]) |
| 246 | + |
| 247 | + cc_obs <- obs_summary$cc[[g]][[r]] |
| 248 | + T_obs <- T_obs + li_cc_stat(cc_obs$r, r_model, cc_obs$n) |
| 249 | + |
| 250 | + ok <- complete.cases(Xg_rep[, c(iv, jv)]) |
| 251 | + r_rep <- cor(Xg_rep[ok, iv], Xg_rep[ok, jv]) |
| 252 | + T_rep <- T_rep + li_cc_stat(r_rep, r_model, sum(ok)) |
| 253 | + } |
| 254 | + } |
| 255 | + } |
| 256 | + |
| 257 | + c(T_obs = T_obs, T_rep = T_rep) |
| 258 | +} |
| 259 | + |
| 260 | +## top-level entry point: limited-information posterior predictive p-value |
| 261 | +## for single-level (single- or multi-group) ordinal/mixed stan/cmdstan |
| 262 | +## fits. Mirrors postpred()'s scaffolding and return shape. |
| 263 | +pp_limited_info <- function(lavobject, thin = 1, parallel = FALSE) { |
| 264 | + lavpartable <- lavobject@ParTable |
| 265 | + lavmodel <- lavobject@Model |
| 266 | + lavoptions <- lavobject@Options |
| 267 | + lavsamplestats <- lavobject@SampleStats |
| 268 | + lavdata <- lavobject@Data |
| 269 | + lavcache <- lavobject@Cache |
| 270 | + lavjags <- lavobject@external$mcmcout |
| 271 | + |
| 272 | + lavmcmc <- make_mcmc(lavjags) |
| 273 | + n.chains <- length(lavmcmc) |
| 274 | + samp.indices <- sampnums(lavjags, thin = thin) |
| 275 | + psamp <- length(samp.indices) |
| 276 | + |
| 277 | + pairs <- li_pairtypes(lavdata, lavsamplestats) |
| 278 | + has_oo <- nrow(pairs$oo) > 0 |
| 279 | + obs_summary <- li_obs_summary(lavdata, pairs) |
| 280 | + |
| 281 | + loop.args <- list(X = 1:n.chains, FUN = function(j) { |
| 282 | + sapply(1:psamp, function(i) { |
| 283 | + li_draw_stat(lavmcmc[[j]][samp.indices[i], ], lavobject, lavmodel, |
| 284 | + lavdata, lavsamplestats, lavpartable, lavoptions, lavcache, |
| 285 | + samp.indices[i], j, pairs, obs_summary, has_oo) |
| 286 | + }) |
| 287 | + }) |
| 288 | + |
| 289 | + if (parallel) { |
| 290 | + loop.args <- c(loop.args, future.seed = TRUE) |
| 291 | + res <- do.call("future_lapply", loop.args) |
| 292 | + } else { |
| 293 | + res <- do.call("lapply", loop.args) |
| 294 | + } |
| 295 | + |
| 296 | + T_obs <- unlist(lapply(res, function(x) x["T_obs", ])) |
| 297 | + T_rep <- unlist(lapply(res, function(x) x["T_rep", ])) |
| 298 | + |
| 299 | + ## draws with an out-of-range implied correlation (see li_oo_stat()) are |
| 300 | + ## NA and excluded from the ppp average, matching how other blavaan |
| 301 | + ## per-draw computations (e.g. get_ll_ord()) drop try()-error draws |
| 302 | + nna <- sum(is.na(T_obs) | is.na(T_rep)) |
| 303 | + if (nna > 0) { |
| 304 | + warning("blavaan WARNING: ", nna, " of ", length(T_obs), |
| 305 | + " posterior draws produced an out-of-range model-implied ", |
| 306 | + "correlation and were excluded from the ppp computation.", |
| 307 | + call. = FALSE) |
| 308 | + } |
| 309 | + |
| 310 | + list(ppval = mean(T_rep > T_obs, na.rm = TRUE), |
| 311 | + ppdist = list(obs = T_obs, reps = T_rep)) |
| 312 | +} |
0 commit comments