--- title: "Heterogeneity: random slopes, latent classes and quantiles" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Heterogeneity: random slopes, latent classes and quantiles} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4) has_qr <- requireNamespace("quantreg", quietly = TRUE) ``` A single coefficient describes the average effect. Health data are often more heterogeneous: patients decline at different rates, a small group of heavy users drives most of the costs, and a treatment can shift the upper tail of a cost distribution more than its centre. This vignette covers three ways to model this heterogeneity: | Heterogeneity | Example | Tool | |---|---|---| | Continuous, in the effect of a variable | each patient has their own rate of lung function decline | random slopes, `panglm(random = ~ time)` | | Discrete, in the whole response | low users and high users of care | latent classes, `panfmm()` | | Across the outcome distribution | effect of chronic disease on high costs | quantile regression, `panqr()` | The examples use small simulated panels and few quadrature nodes, EM starts and bootstrap replicates so that the vignette builds quickly; the companion book uses larger settings. ```{r setup} library(panglm) ``` ## Random slopes A random-intercept model lets patients differ in their level only. A random-slope model also lets them differ in the effect of a regressor, usually time: $y_{it} = x_{it}^{\mathsf T}\beta + b_{0i} + b_{1i} t + \varepsilon_{it}$, with $(b_{0i}, b_{1i})$ bivariate normal with an unstructured covariance (Laird and Ware 1982). We simulate lung function (FEV1) of 100 patients over 6 years, declining by 50 mL a year on average, 30 mL faster in smokers, with a patient-specific rate of decline (SD 40 mL a year). ```{r rs-data} set.seed(3) N <- 100; TT <- 6 f <- data.frame(id = rep(seq_len(N), each = TT), year = rep(0:(TT - 1), N)) f$smoker <- rep(rbinom(N, 1, 0.4), each = TT) b0 <- rnorm(N, sd = 0.4) b1 <- rnorm(N, sd = 0.04) f$fev1 <- 3 + b0[f$id] + (-0.05 - 0.03 * f$smoker + b1[f$id]) * f$year - 0.2 * f$smoker + rnorm(N * TT, sd = 0.1) f$exac <- rpois(N * TT, exp(-0.3 + 0.3 * f$smoker + (0.05 + 3 * b1[f$id]) * f$year + b0[f$id])) ``` `random = ~ year` adds a random slope of `year` to the random intercept. The Gaussian model is fitted by exact maximum likelihood. ```{r rs} rs <- panglm(fev1 ~ year * smoker, data = f, index = c("id", "year"), model = "random", random = ~ year) summary(rs) head(rs$ranef, 3) ``` The average decline is 47 mL a year, 17 mL faster in smokers. The SD of the patient-specific slopes, 42 mL a year, is as large as the average decline itself: some patients barely decline, others lose about 130 mL a year. `rs$varcomp` holds the covariance matrix of the random effects and `rs$ranef` the predicted effects of each patient (BLUPs). A likelihood ratio test against the random-intercept model rejects clearly. With a variance on the boundary of its space under the null, the $\chi^2_2$ p-value is conservative (Self and Liang 1987). ```{r rs-lrt} ri <- panglm(fev1 ~ year * smoker, data = f, index = c("id", "year"), model = "random") lr <- 2 * (as.numeric(logLik(rs)) - as.numeric(logLik(ri))) c(LR = lr, p = pchisq(lr, df = 2, lower.tail = FALSE)) ``` Random slopes are also available for Poisson, negative binomial, gamma and binomial models. Their likelihood is integrated by adaptive Gauss-Hermite quadrature with `R` nodes per random effect (`R = 1` is the Laplace approximation of `lme4::glmer()`). Use at least 5 nodes and check that the estimates are stable when `R` is increased. ```{r rs-pois} rp <- panglm(exac ~ year + smoker, data = f, index = c("id", "year"), model = "random", family = "poisson", random = ~ year, R = 5) summary(rp) ``` The exacerbation rate rises by about 10% a year on average; the slope SD of 0.12 means that for some patients it falls and for others it rises by more than 30% a year. ## Latent classes Finite mixture models assume that each individual belongs to one of $k$ latent classes for all of their observations, with class-specific coefficients (Deb and Trivedi 1997, 2002). In health care use the classes typically separate a large group of low users from a small group of high users who respond differently to illness or insurance. `concomitant` lets the class probabilities depend on individual-level variables through a multinomial logit. We simulate 300 patients, 30% of them high users with a stronger effect of a chronic condition. High users are more frequent among patients with low income. ```{r fmm-data} set.seed(1) N <- 300; TT <- 4 m <- data.frame(id = rep(seq_len(N), each = TT), year = rep(seq_len(TT), N)) income <- rnorm(N) high <- rbinom(N, 1, plogis(-1 - 0.8 * income)) m$income <- income[m$id] m$chronic <- rbinom(N * TT, 1, 0.4) m$age <- rep(round(runif(N, 40, 80)), each = TT) m$visits <- rpois(N * TT, exp(ifelse(high[m$id] == 1, 1.5 + 0.6 * m$chronic, 0 + 0.3 * m$chronic) + 0.01 * (m$age - 60))) ``` The number of classes is chosen first, usually by BIC (Deb and Trivedi 2002 prefer it). Mixture likelihoods have local maxima, so the EM algorithm is run from several starts (`nstart`). The comparison below uses a single start and the final model three, to keep the example fast; in practice use 10 or more. EM converges slowly when a model has more classes than the data support; a looser tolerance (`tol = 1e-4`) is enough to compare the fits and leaves this table unchanged. ```{r fmm-select} panfmmselect(visits ~ chronic + age, data = m, index = c("id", "year"), k = 1:3, nstart = 1, seed = 1, tol = 1e-4) ``` BIC picks two classes; the third class of the three-class model holds only 5% of the patients. The final model adds income as a concomitant variable. ```{r fmm} fmm <- panfmm(visits ~ chronic + age, data = m, index = c("id", "year"), k = 2, family = "poisson", concomitant = ~ income, nstart = 3, seed = 1) summary(fmm) ``` The model finds a class of about 70% low users and a class of high users with more than five times as many visits and a larger effect of chronic illness (a rate ratio of $e^{0.65} = 1.9$ against $e^{0.31} = 1.4$; true values 1.8 and 1.35). Higher income lowers the probability of being a high user (true coefficient -0.8). With classes this well separated, every patient is assigned to their true class. `predict(type = "posterior")` gives each patient's posterior class probabilities; `type = "class"` assigns each to the class with the largest posterior probability, and `type = "components"` gives the class-specific means. ```{r fmm-post, message = FALSE} first <- !duplicated(m$id) post <- predict(fmm, type = "posterior")[first, ] assigned <- predict(fmm, type = "class")[first] all(assigned == max.col(post)) table(true = high, assigned = assigned) ``` Latent classes are a modelling device, not a discovery of real groups: a two-class mixture also fits a continuous but skewed heterogeneity. Interpret the classes with care and check that the conclusions hold with $k + 1$ classes. ## Quantile regression Mean regression can hide what happens in the tails. Quantile regression estimates the effect of a regressor on each quantile $\tau$ of the outcome distribution (Koenker and Bassett 1978). With panel data the individual effects have to be removed, which is harder than for the mean because quantiles do not commute with differencing. `panqr()` offers the two-step estimator of Canay (2011), which subtracts individual effects estimated at the mean, and the method of moments quantile regression (MM-QR) of Machado and Santos Silva (2019), which assumes a location-scale model and is the estimator of Stata's `xtqreg`. Both assume that the individual effects shift the whole distribution (a location shift). We simulate annual costs whose spread grows with the number of chronic conditions, so that the effect of chronic conditions is larger in the upper tail. ```{r qr-data} set.seed(1) N <- 200; TT <- 5 q <- data.frame(id = rep(seq_len(N), each = TT), year = rep(seq_len(TT), N)) alpha <- rep(rnorm(N), each = TT) q$chronic <- rpois(N * TT, 1 + 0.5 * (alpha > 0)) q$age <- rep(sample(40:70, N, TRUE), each = TT) + q$year q$cost <- 2 + alpha + 0.5 * q$chronic + 0.02 * q$age + (1 + 0.3 * q$chronic) * rexp(N * TT) ``` ```{r qr, eval = has_qr} qr <- panqr(cost ~ chronic + age, data = q, index = c("id", "year"), tau = c(0.1, 0.25, 0.5, 0.75, 0.9), model = "mmqr") summary(qr) ``` ```{r qr-plot, eval = has_qr, fig.alt = "Coefficient of chronic and age across quantiles with confidence bands"} plot(qr) ``` Each additional chronic condition raises the median cost by about 0.7 and the 90th percentile by about 0.9 (true values 0.71 and 1.19): the effect grows across the distribution, as the significant scale coefficient indicates. The age effect is small and imprecise at every quantile. The MM-QR standard errors are clustered by individual, so they allow the costs of a patient to be correlated across years. The warning flags the few observations at which the fitted scale is not positive, where the location-scale model is strained. The Canay estimator uses a cluster bootstrap for its standard errors (50 replicates here, use 500 or more); a pooled quantile regression that ignores the individual effects is biased when, as here, they are correlated with the regressors (0.85 and 1.23 below). ```{r qr-canay, eval = has_qr, warning = FALSE} cn <- panqr(cost ~ chronic + age, data = q, index = c("id", "year"), tau = c(0.5, 0.9), model = "canay", boot = 50, seed = 1) coef(cn) coef(panqr(cost ~ chronic + age, data = q, index = c("id", "year"), tau = c(0.5, 0.9), model = "pooling")) ``` ## References Canay, I. A. (2011). A simple approach to quantile regression for panel data. *The Econometrics Journal*, 14, 368-386. Deb, P., and Trivedi, P. K. (1997). Demand for medical care by the elderly: a finite mixture approach. *Journal of Applied Econometrics*, 12, 313-336. Deb, P., and Trivedi, P. K. (2002). The structure of demand for health care: latent class versus two-part models. *Journal of Health Economics*, 21, 601-625. Koenker, R., and Bassett, G. (1978). Regression quantiles. *Econometrica*, 46, 33-50. Laird, N. M., and Ware, J. H. (1982). Random-effects models for longitudinal data. *Biometrics*, 38, 963-974. Machado, J. A. F., and Santos Silva, J. M. C. (2019). Quantiles via moments. *Journal of Econometrics*, 213, 145-173. Self, S. G., and Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. *Journal of the American Statistical Association*, 82, 605-610.