Processing math: 100%

Longitudinal hurdle models 3

In the last post on longitudinal hurdle models, I had just taken samples from the marginal mean

g(θ;x)=E(Yθ;x)=E(Yui,vi,θ;x)f(ui,viθ,y)duidviL1Ll=1E(Yu(l)i,v(l)i,θ;x)=L1Ll=1{1logit1(xTγ+u(l)i)}exp(xTβ+v(l)i+σ22).

One of the issues with lognormal data is that it is highly skewed, so the mean can be very large. In a small sample, the sample mean can change a lot based on just 1 or 2 large observations. For this reason I would like to sample from other summary measures of Y.

Samples from p(Ykθ;x)

This is quite similar to taking samples from the marginal mean.

r(θ;x,k)=p(Y<kθ;x)=p(Y<kui,vi,θ;x)f(ui,viθ,y)duidviL1Ll=1p(Y<ku(l)i,v(l)i,θ;x)=L1Ll=1 [π(l)(x)+{1π(l)(x)}Φ{log(k)xTβv(l)iσ}].

where π(l)(x):=logit1(xTγ+u(l)i).

Again, I don’t know of any functions for doing this, so I built my own.

devtools::install_github("dominicmagirr/hurlong")
library(brms)
x <- data.frame(id = NA, time = "2")
hurlong::marg_pyk_q(k = 0.5,
                    newdata = x, 
                    nsims = 1000, 
                    fit = fit_hurdle)
##   id time     2.5%       50%     97.5%
## 1 NA    2 0.421836 0.5018633 0.5792536

Samples from median(Yθ;x)

Finally, I am interested in samples from the marginal median.

s(θ;x)=median(Yθ;x)=median(Yui,vi,θ;x)f(ui,viθ,y)duidviL1Ll=1median(Yu(l)i,v(l)i,θ;x)

To evaluate median(Yu(l)i,v(l)i,θ;x) at each l, I do a search for m(l) such that

p(Y<m(l)u(l)i,v(l)i,θ;x)=0.5.

Again, I’ve made a function that can do this (for this specific longitudinal hurdle model).

hurlong::marg_med(newdata = x, 
                  nsims = 1000, 
                  ks = exp(seq(log(0.01), 
                               log(100), 
                               length.out = 15)), # where to evaluate p(Y<k)
                  fit = fit_hurdle)
##   id time      2.5%       50%     97.5%
## 1 NA    2 0.2332964 0.4894882 0.8856152

Predictions

To round off this series on longitudinal hurdle models, I want to show how to simulate draws (and find quantiles) from the posterior predictive distribution for a new observation (˜Y). Firstly for a patient i already in the data set, where we draw from

f(˜yy;x)=f(˜yui,vi,θ,y;x)f(ui,vi,θy)

library(brms)
predict(fit_hurdle,
        newdata = data.frame(id = 1, time = "2"),
        robust = TRUE)
##       Estimate Est.Error Q2.5    Q97.5
## [1,] 0.2448615 0.3224153    0 2.145083

and, secondly, for a new patient i, where we draw from

f(˜yy;x)=f(˜yui,vi,θ,y;x)f(ui,viθ,y)f(θy)

predict(fit_hurdle,
        newdata = data.frame(id = NA, time = "2"),
        allow_new_levels = TRUE,
        sample_new_levels = "gaussian",
        robust = TRUE)
##       Estimate Est.Error Q2.5    Q97.5
## [1,] 0.4614108 0.6840876    0 52.57643
Avatar
Dominic Magirr
Medical Statistician

Interested in the design, analysis and interpretation of clinical trials.

Related