Data set-up and calculation of marginal probabilities from a continuation ratio model

cr_setup(y, direction = c("forward", "backward"))

cr_marg_probs(eta, direction = c("forward", "backward"))

## Arguments

y

a numeric vector denoting the ordinal response variable.

direction

character string specifying the direction of the continuation ratio model; "forward" corresponds to a discrete hazard function.

eta

a numeric matrix of the linear predictor, with columns corresponding to the different levels of the ordinal response.

## Author

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

Frank Harrell

## Note

Function cr_setup() is based on the cr.setup() function from package rms.

## Examples

n <- 300 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 15 # maximum follow-up time

# we constuct a data frame with the design:
# everyone has a baseline measurment, and then measurements at random follow-up times
DF <- data.frame(id = rep(seq_len(n), each = K),
time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))

# design matrices for the fixed and random effects
X <- model.matrix(~ sex * time, data = DF)[, -1]
Z <- model.matrix(~ 1, data = DF)

thrs <- c(-1.5, 0, 0.9) # thresholds for the different ordinal categories
betas <- c(-0.25, 0.24, -0.05) # fixed effects coefficients
D11 <- 0.48 # variance of random intercepts
D22 <- 0.1 # variance of random slopes

# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))[, 1, drop = FALSE]
# linear predictor
eta_y <- drop(X %*% betas + rowSums(Z * b[DF$id, , drop = FALSE])) # linear predictor for each category eta_y <- outer(eta_y, thrs, "+") # marginal probabilities per category mprobs <- cr_marg_probs(eta_y) # we simulate ordinal longitudinal data DF$y <- unname(apply(mprobs, 1, sample, x = ncol(mprobs), size = 1, replace = TRUE))

# If you want to simulate from the backward formulation of the CR model, you need to
# change eta_y <- outer(eta_y, thrs, "+") to eta_y <- outer(eta_y, rev(thrs), "+"),
# and mprobs <- cr_marg_probs(eta_y) to mprobs <- cr_marg_probs(eta_y, "backward")

#################################################

# prepare the data
# If you want to fit the CR model under the backward formulation, you need to change
# cr_vals <- cr_setup(DF$y) to cr_vals <- cr_setup(DF$y, "backward")
cr_vals <- cr_setup(DF$y) cr_data <- DF[cr_vals$subs, ]
cr_data$y_new <- cr_vals$y
cr_data$cohort <- cr_vals$cohort

# fit the model
fm <- mixed_model(y_new ~ cohort + sex * time, random = ~ 1 | id,
data = cr_data, family = binomial())

summary(fm)
#>
#> Call:
#> mixed_model(fixed = y_new ~ cohort + sex * time, random = ~1 |
#>     id, data = cr_data, family = binomial())
#>
#> Data Descriptives:
#> Number of Observations: 4216
#> Number of Groups: 300
#>
#> Model:
#>  family: binomial
#>
#> Fit statistics:
#>    log.Lik      AIC      BIC
#>  -2401.604 4817.208 4843.134
#>
#> Random effects covariance matrix:
#>               StdDev
#> (Intercept) 0.673666
#>
#> Fixed effects:
#>                Estimate Std.Err  z-value   p-value
#> (Intercept)     -1.6576  0.1138 -14.5648   < 1e-04
#> cohorty>=2       1.6832  0.0909  18.5238   < 1e-04
#> cohorty>=3       2.5602  0.1390  18.4215   < 1e-04
#> sexfemale       -0.3591  0.1410  -2.5467 0.0108744
#> time             0.2590  0.0139  18.6439   < 1e-04
#> sexfemale:time  -0.0473  0.0172  -2.7413 0.0061192
#>
#> Integration: