rvlenth / emmeans

Estimated marginal means
https://rvlenth.github.io/emmeans/
340 stars 30 forks source link

Lost trend (ordered factor) when using qdrg in the multiple imputation context. #487

Closed Generalized closed 2 months ago

Generalized commented 2 months ago

Background: I have a longitudinal study. I want to assess trend in LS-means. Variable "visit" is ordinal factor with Levels: Month 6 < Month 12 < Month 20 I also want to use other than default estimator of covariance.

Problem: emmeans(as.mira(...)) takes defaults and runs fine. But I want to provide my own settings. So I need the qdrg() to use pooled myself parameters and vcov (so I can use the one I need). But the qrdg() way misses some things and it doesn't work.

Examples:

Let's fit the model first:

> imp_mod <- lapply(d, 
                  function(data)
                    glmgee(ODIPain ~ Visit * Arm,
                    id = PatientId,
                      data = data,
                      family = gaussian(link = "identity"),
                      corstr = "unstructured")) 

> imp_mod 
$`1`

        Variance function:  gaussian
                     Link:  identity
    Correlation structure:  Unstructured 

$`2`

        Variance function:  gaussian
                     Link:  identity
    Correlation structure:  Unstructured 

> summary(imp_mod[[1]])

Sample size
   Number of observations:  333
       Number of clusters:  111 
             Cluster size:  3 
*************************************************************
Model
        Variance function:  gaussian
            Link function:  identity
    Correlation structure:  Unstructured
*************************************************************
Coefficients
             Estimate Std.Error z-value Pr(>|z|)
(Intercept)     1.903     0.158  12.029   <2e-16
Visit.L         0.244     0.139   1.761   0.0783
Visit.Q        -0.096     0.130  -0.741   0.4588
ArmB           -0.230     0.208  -1.109   0.2676
Visit.L:ArmB   -0.459     0.187  -2.459   0.0140
Visit.Q:ArmB    0.045     0.167   0.272   0.7856

Dispersion      1.807                           
*************************************************************
Working correlation
     [1]   [2]   [3] 
[1] 1.000 0.472 0.428
[2] 0.472 1.000 0.691
[3] 0.428 0.691 1.000

So far the "L" and "Q" is recognized by the fitting procedure.

Now let's pool and use basic emmeans over as.mira(...):

> (imp_mod_em <- emmeans(as.mira(imp_mod), 
+                        specs = ~ Visit | Arm, adjust="none"))
Arm = A:
 Visit    emmean    SE  df lower.CL upper.CL
 Month 6    1.73 0.184 327     1.36     2.09
 Month 12   1.96 0.207 327     1.55     2.36
 Month 20   2.06 0.203 327     1.66     2.45

Arm = B:
 Visit    emmean    SE  df lower.CL upper.CL
 Month 6    1.80 0.165 327     1.48     2.13
 Month 12   1.75 0.187 327     1.38     2.12
 Month 20   1.54 0.171 327     1.20     1.87

Covariance estimate used: robust 
Confidence level used: 0.95 

and test the trends:

> update(contrast(imp_mod_em, "poly", by = "Arm"), 
+        adjust="none", level = 0.95, infer = c(TRUE, TRUE))
Arm = A:
 contrast  estimate    SE  df lower.CL upper.CL t.ratio p.value
 linear       0.327 0.204 327   -0.074    0.729   1.603  0.1100
 quadratic   -0.127 0.376 327   -0.867    0.612  -0.339  0.7350

Arm = B:
 contrast  estimate    SE  df lower.CL upper.CL t.ratio p.value
 linear      -0.268 0.182 327   -0.626    0.090  -1.472  0.1420
 quadratic   -0.161 0.272 327   -0.696    0.375  -0.590  0.5550

Confidence level used: 0.95 

OK, it works well.


Now, I don't want the default sandwich estimator of covariance but rather the bias-corrected robust one. I need to pass somehow this parameter, so I need to use qdrg and some programming.

This function will be required. It adapts the code from emmeans to pool estimates:

pool_estimates_for_qdrg <- function(k, coefficients, varcov) {

  V <- 1/k * varcov[[1]]
  allb <- cbind(coefficients[[1]], matrix(0, nrow = length(coefficients[[1]]), ncol = k - 1))

  for (i in 1 + seq_len(k - 1)) {
    V <- V + 1/k * varcov[[i]]
    allb[, i] <- coefficients[[i]]
  }

  pooled_coef <- apply(allb, 1, mean)
  notna <- which(!is.na(pooled_coef))
  pooled_VC <- V + (k + 1)/k * cov(t(allb[notna, , drop = FALSE])) 

  return(list(pooled_coef = pooled_coef, pooled_VC = pooled_VC))
}

Let's pool the coefficients:

> (pooled_coefficients <- pool_estimates_for_qdrg(k=2, 
+                         coefficients = lapply(imp_mod, function(mod) mod$coefficients),
+                         varcov       = lapply(imp_mod, function(mod) vcov(mod, type = "bias-corrected"))))
$pooled_coef
 (Intercept)      Visit.L      Visit.Q         ArmB Visit.L:ArmB Visit.Q:ArmB 
      1.9121       0.2314      -0.0520      -0.2157      -0.4208      -0.0137 

$pooled_VC
             (Intercept)  Visit.L  Visit.Q     ArmB Visit.L:ArmB Visit.Q:ArmB
(Intercept)      0.02542  0.00363 -0.00161 -0.02477     -0.00294      0.00121
Visit.L          0.00363  0.02160 -0.00227 -0.00455     -0.02258      0.00283
Visit.Q         -0.00161 -0.00227  0.02421  0.00479      0.00565     -0.02616
ArmB            -0.02477 -0.00455  0.00479  0.04559      0.00591     -0.00809
Visit.L:ArmB    -0.00294 -0.02258  0.00565  0.00591      0.04066     -0.00935
Visit.Q:ArmB     0.00121  0.00283 -0.02616 -0.00809     -0.00935      0.04089

Now the grid:

> imp_mod_grid <- with(pooled_coefficients,
+                      qdrg(formula = ODIPain ~ Visit * Arm,
+                           data=d[[1]], # fake data in any case 
+                           coef = pooled_coef, 
+                           vcov = pooled_VC))
> (imp_mod_grid <- with(pooled_coefficients,
+                      qdrg(formula = ODIPain ~ Visit * Arm,
+                           data=d[[1]], # fake data in any case 
+                           coef = pooled_coef, 
+                           vcov = pooled_VC)))
'emmGrid' object with variables:
    Visit = Month 6, Month 12, Month 20
    Arm = A, B

When passed to emmeans, all values are same:

(imp_mod_em <- emmeans(imp_mod_grid, 
                       specs = ~ Visit | Arm, adjust="none"))

Arm = A:
 Visit    emmean    SE  df asymp.LCL asymp.UCL
 Month 6    1.91 0.159 Inf      1.60      2.23
 Month 12   1.91 0.159 Inf      1.60      2.23
 Month 20   1.91 0.159 Inf      1.60      2.23

Arm = B:
 Visit    emmean    SE  df asymp.LCL asymp.UCL
 Month 6    1.70 0.146 Inf      1.41      1.98
 Month 12   1.70 0.146 Inf      1.41      1.98
 Month 20   1.70 0.146 Inf      1.41      1.98

Confidence level used: 0.95 

and the trends are "empty":

> update(contrast(imp_mod_em, "poly", by = "Arm"), 
+        adjust="none", level = 0.95, infer = c(TRUE, TRUE))
Arm = A:
 contrast  estimate SE  df asymp.LCL asymp.UCL z.ratio p.value
 linear           0  0 Inf         0         0     NaN     NaN
 quadratic        0  0 Inf         0         0     NaN     NaN

Arm = B:
 contrast  estimate SE  df asymp.LCL asymp.UCL z.ratio p.value
 linear           0  0 Inf         0         0     NaN     NaN
 quadratic        0  0 Inf         0         0     NaN     NaN

Confidence level used: 0.95 

Is it possible to obtain the estimates for trends in this setting? I need the qdrg() to use specific covariance estimator and there seems to be no other way where I can use it. But somehow the information about the ordering in visits is lost or something else happens...


The data: 2 imputed datasets in a list.

d <- list(`1` = structure(list(PatientId = c("ccbbaa04", "ccbbaa04", 
"ccbbaa04", "ccbbaa07", "ccbbaa07", "ccbbaa07", "ccbbaa08", "ccbbaa08", 
"ccbbaa08", "ccbbaa10", "ccbbaa10", "ccbbaa10", "ccbbaa12", "ccbbaa12", 
"ccbbaa12", "ccbbaa14", "ccbbaa14", "ccbbaa14", "ccbbaa17", "ccbbaa17", 
"ccbbaa17", "ccbbaa18", "ccbbaa18", "ccbbaa18", "ccbbaa19", "ccbbaa19", 
"ccbbaa19", "ccbbaa21", "ccbbaa21", "ccbbaa21", "ccbbaa23", "ccbbaa23", 
"ccbbaa23", "ccbbaa24", "ccbbaa24", "ccbbaa24", "ccbbaa29", "ccbbaa29", 
"ccbbaa29", "ccbbaa30", "ccbbaa30", "ccbbaa30", "ccbbaa31", "ccbbaa31", 
"ccbbaa31", "ccbbaa33", "ccbbaa33", "ccbbaa33", "ccbbaa34", "ccbbaa34", 
"ccbbaa34", "ccbbaa36", "ccbbaa36", "ccbbaa36", "ccbbaa37", "ccbbaa37", 
"ccbbaa37", "ccbbaa38", "ccbbaa38", "ccbbaa38", "ccbbaa41", "ccbbaa41", 
"ccbbaa41", "ccbbaa42", "ccbbaa42", "ccbbaa42", "ccbbaa43", "ccbbaa43", 
"ccbbaa43", "ccbbaa44", "ccbbaa44", "ccbbaa44", "ccbbaa45", "ccbbaa45", 
"ccbbaa45", "ccbbaa47", "ccbbaa47", "ccbbaa47", "ccbbaa48", "ccbbaa48", 
"ccbbaa48", "ccbbaa49", "ccbbaa49", "ccbbaa49", "ccbbaa50", "ccbbaa50", 
"ccbbaa50", "ccbbaa51", "ccbbaa51", "ccbbaa51", "ccbbaa52", "ccbbaa52", 
"ccbbaa52", "ccbbaa54", "ccbbaa54", "ccbbaa54", "ccbbaa55", "ccbbaa55", 
"ccbbaa55", "ccbbaa56", "ccbbaa56", "ccbbaa56", "ccbbaa57", "ccbbaa57", 
"ccbbaa57", "ccbbaa58", "ccbbaa58", "ccbbaa58", "ccbbaa59", "ccbbaa59", 
"ccbbaa59", "ccbbaa61", "ccbbaa61", "ccbbaa61", "ccbbaa62", "ccbbaa62", 
"ccbbaa62", "ccbbaa63", "ccbbaa63", "ccbbaa63", "ccbbaa64", "ccbbaa64", 
"ccbbaa64", "ccbbaa66", "ccbbaa66", "ccbbaa66", "ccbbaa68", "ccbbaa68", 
"ccbbaa68", "ccbbaa69", "ccbbaa69", "ccbbaa69", "ccbbaa70", "ccbbaa70", 
"ccbbaa70", "ccbbaa71", "ccbbaa71", "ccbbaa71", "ccbbaa72", "ccbbaa72", 
"ccbbaa72", "ccbbaa73", "ccbbaa73", "ccbbaa73", "ccbbaa74", "ccbbaa74", 
"ccbbaa74", "ccbbaa75", "ccbbaa75", "ccbbaa75", "ccbbaa76", "ccbbaa76", 
"ccbbaa76", "ccbbaa77", "ccbbaa77", "ccbbaa77", "ccddaabb", "ccddaabb", 
"ccddaabb", "ccddaacc", "ccddaacc", "ccddaacc", "ccddaa05", "ccddaa05", 
"ccddaa05", "ccddaa09", "ccddaa09", "ccddaa09", "ccddaa14", "ccddaa14", 
"ccddaa14", "ccddaa15", "ccddaa15", "ccddaa15", "ccddaa17", "ccddaa17", 
"ccddaa17", "ccddaa20", "ccddaa20", "ccddaa20", "ccddaa21", "ccddaa21", 
"ccddaa21", "ccddaa22", "ccddaa22", "ccddaa22", "ccddaa26", "ccddaa26", 
"ccddaa26", "ccddaa27", "ccddaa27", "ccddaa27", "ccddaa28", "ccddaa28", 
"ccddaa28", "ccddaa29", "ccddaa29", "ccddaa29", "ccddaa30", "ccddaa30", 
"ccddaa30", "ccddaa31", "ccddaa31", "ccddaa31", "ccddaa32", "ccddaa32", 
"ccddaa32", "ccddaa33", "ccddaa33", "ccddaa33", "ccddaa34", "ccddaa34", 
"ccddaa34", "ccddaa35", "ccddaa35", "ccddaa35", "ccddaa36", "ccddaa36", 
"ccddaa36", "ccddaa37", "ccddaa37", "ccddaa37", "ccddaa38", "ccddaa38", 
"ccddaa38", "ccddaa39", "ccddaa39", "ccddaa39", "ccddaa40", "ccddaa40", 
"ccddaa40", "ccddaa41", "ccddaa41", "ccddaa41", "ccddaa42", "ccddaa42", 
"ccddaa42", "ccddaa43", "ccddaa43", "ccddaa43", "ccddaa45", "ccddaa45", 
"ccddaa45", "ccddaa47", "ccddaa47", "ccddaa47", "ccddaa48", "ccddaa48", 
"ccddaa48", "ccddaa49", "ccddaa49", "ccddaa49", "ccddaa50", "ccddaa50", 
"ccddaa50", "ccddaa51", "ccddaa51", "ccddaa51", "ccddaa52", "ccddaa52", 
"ccddaa52", "ccccaabb", "ccccaabb", "ccccaabb", "ccccaa04", "ccccaa04", 
"ccccaa04", "ccccaa05", "ccccaa05", "ccccaa05", "ccccaa06", "ccccaa06", 
"ccccaa06", "ccccaa09", "ccccaa09", "ccccaa09", "ccccaa12", "ccccaa12", 
"ccccaa12", "ccccaa13", "ccccaa13", "ccccaa13", "ccccaa14", "ccccaa14", 
"ccccaa14", "ccccaa15", "ccccaa15", "ccccaa15", "ccccaa16", "ccccaa16", 
"ccccaa16", "ccccaa17", "ccccaa17", "ccccaa17", "ccccaa18", "ccccaa18", 
"ccccaa18", "ccccaa19", "ccccaa19", "ccccaa19", "ccccaa20", "ccccaa20", 
"ccccaa20", "ccccaa22", "ccccaa22", "ccccaa22", "ccccaa23", "ccccaa23", 
"ccccaa23", "ccccaa24", "ccccaa24", "ccccaa24", "ccccaa25", "ccccaa25", 
"ccccaa25", "ccccaa26", "ccccaa26", "ccccaa26", "ccccaa27", "ccccaa27", 
"ccccaa27", "ccccaa28", "ccccaa28", "ccccaa28", "ccccaa30", "ccccaa30", 
"ccccaa30", "ccccaa31", "ccccaa31", "ccccaa31", "ccccaa32", "ccccaa32", 
"ccccaa32"), Arm = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 
1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 
1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 
2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 
1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L), levels = c("A", "B", "C"), class = "factor"), 
    Visit = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), levels = c("Month 6", "Month 12", 
    "Month 20"), class = c("ordered", "factor")), ODIPain = c(2L, 
    1L, 1L, 0L, 0L, 2L, 1L, 3L, 0L, 2L, 0L, 0L, 3L, 1L, 1L, 2L, 
    1L, 2L, 3L, 3L, 3L, 1L, 2L, 1L, 2L, 2L, 3L, 4L, 2L, 1L, 1L, 
    0L, 0L, 4L, 3L, 3L, 4L, 2L, 1L, 5L, 0L, 1L, 1L, 2L, 2L, 0L, 
    2L, 3L, 0L, 2L, 3L, 1L, 1L, 1L, 2L, 1L, 3L, 2L, 2L, 1L, 0L, 
    0L, 0L, 2L, 5L, 4L, 2L, 2L, 3L, 2L, 4L, 5L, 1L, 4L, 3L, 4L, 
    4L, 4L, 4L, 1L, 3L, 1L, 2L, 3L, 4L, 2L, 2L, 2L, 2L, 2L, 4L, 
    2L, 2L, 0L, 3L, 3L, 2L, 2L, 3L, 2L, 4L, 0L, 2L, 2L, 2L, 0L, 
    0L, 1L, 4L, 3L, 4L, 2L, 3L, 3L, 4L, 3L, 3L, 0L, 1L, 2L, 0L, 
    1L, 0L, 1L, 4L, 4L, 0L, 0L, 0L, 1L, 0L, 0L, 2L, 2L, 2L, 1L, 
    1L, 2L, 2L, 2L, 5L, 0L, 2L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 
    2L, 0L, 3L, 0L, 0L, 5L, 5L, 5L, 3L, 3L, 2L, 2L, 4L, 2L, 2L, 
    1L, 0L, 2L, 0L, 0L, 2L, 0L, 0L, 0L, 3L, 1L, 3L, 4L, 5L, 4L, 
    3L, 2L, 3L, 4L, 4L, 2L, 2L, 0L, 2L, 3L, 2L, 0L, 0L, 2L, 2L, 
    1L, 1L, 0L, 0L, 1L, 3L, 5L, 3L, 2L, 2L, 1L, 2L, 2L, 3L, 1L, 
    2L, 2L, 2L, 3L, 3L, 2L, 2L, 4L, 2L, 2L, 3L, 2L, 0L, 2L, 2L, 
    2L, 2L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 3L, 1L, 0L, 1L, 1L, 
    1L, 1L, 0L, 2L, 0L, 2L, 2L, 2L, 3L, 4L, 2L, 3L, 2L, 2L, 2L, 
    2L, 0L, 1L, 0L, 2L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 
    1L, 2L, 0L, 2L, 0L, 0L, 0L, 0L, 2L, 1L, 1L, 0L, 0L, 0L, 2L, 
    2L, 1L, 0L, 0L, 0L, 2L, 3L, 2L, 3L, 5L, 5L, 3L, 3L, 2L, 2L, 
    4L, 2L, 2L, 2L, 2L, 0L, 0L, 0L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 
    4L, 3L, 2L, 4L, 2L, 2L, 2L, 2L, 4L, 4L, 3L, 1L, 2L, 1L, 1L, 
    0L, 1L)), row.names = c(56L, 167L, 278L, 57L, 168L, 279L, 
1L, 112L, 223L, 58L, 169L, 280L, 2L, 113L, 224L, 59L, 170L, 281L, 
3L, 114L, 225L, 60L, 171L, 282L, 61L, 172L, 283L, 62L, 173L, 
284L, 4L, 115L, 226L, 5L, 116L, 227L, 6L, 117L, 228L, 63L, 174L, 
285L, 64L, 175L, 286L, 7L, 118L, 229L, 8L, 119L, 230L, 65L, 176L, 
287L, 9L, 120L, 231L, 66L, 177L, 288L, 67L, 178L, 289L, 10L, 
121L, 232L, 68L, 179L, 290L, 11L, 122L, 233L, 12L, 123L, 234L, 
69L, 180L, 291L, 13L, 124L, 235L, 14L, 125L, 236L, 70L, 181L, 
292L, 71L, 182L, 293L, 72L, 183L, 294L, 15L, 126L, 237L, 16L, 
127L, 238L, 17L, 128L, 239L, 73L, 184L, 295L, 18L, 129L, 240L, 
19L, 130L, 241L, 74L, 185L, 296L, 20L, 131L, 242L, 75L, 186L, 
297L, 76L, 187L, 298L, 77L, 188L, 299L, 21L, 132L, 243L, 22L, 
133L, 244L, 78L, 189L, 300L, 23L, 134L, 245L, 24L, 135L, 246L, 
79L, 190L, 301L, 25L, 136L, 247L, 80L, 191L, 302L, 26L, 137L, 
248L, 81L, 192L, 303L, 27L, 138L, 249L, 82L, 193L, 304L, 28L, 
139L, 250L, 83L, 194L, 305L, 29L, 140L, 251L, 84L, 195L, 306L, 
85L, 196L, 307L, 86L, 197L, 308L, 87L, 198L, 309L, 30L, 141L, 
252L, 88L, 199L, 310L, 89L, 200L, 311L, 31L, 142L, 253L, 32L, 
143L, 254L, 33L, 144L, 255L, 90L, 201L, 312L, 91L, 202L, 313L, 
92L, 203L, 314L, 34L, 145L, 256L, 35L, 146L, 257L, 36L, 147L, 
258L, 37L, 148L, 259L, 38L, 149L, 260L, 93L, 204L, 315L, 94L, 
205L, 316L, 95L, 206L, 317L, 39L, 150L, 261L, 96L, 207L, 318L, 
97L, 208L, 319L, 40L, 151L, 262L, 98L, 209L, 320L, 41L, 152L, 
263L, 42L, 153L, 264L, 99L, 210L, 321L, 43L, 154L, 265L, 100L, 
211L, 322L, 44L, 155L, 266L, 45L, 156L, 267L, 101L, 212L, 323L, 
46L, 157L, 268L, 102L, 213L, 324L, 47L, 158L, 269L, 48L, 159L, 
270L, 103L, 214L, 325L, 49L, 160L, 271L, 104L, 215L, 326L, 50L, 
161L, 272L, 105L, 216L, 327L, 51L, 162L, 273L, 52L, 163L, 274L, 
106L, 217L, 328L, 107L, 218L, 329L, 53L, 164L, 275L, 54L, 165L, 
276L, 108L, 219L, 330L, 55L, 166L, 277L, 109L, 220L, 331L, 110L, 
221L, 332L, 111L, 222L, 333L), class = "data.frame"), `2` = structure(list(
    PatientId = c("ccbbaa04", "ccbbaa04", "ccbbaa04", "ccbbaa07", 
    "ccbbaa07", "ccbbaa07", "ccbbaa08", "ccbbaa08", "ccbbaa08", 
    "ccbbaa10", "ccbbaa10", "ccbbaa10", "ccbbaa12", "ccbbaa12", 
    "ccbbaa12", "ccbbaa14", "ccbbaa14", "ccbbaa14", "ccbbaa17", 
    "ccbbaa17", "ccbbaa17", "ccbbaa18", "ccbbaa18", "ccbbaa18", 
    "ccbbaa19", "ccbbaa19", "ccbbaa19", "ccbbaa21", "ccbbaa21", 
    "ccbbaa21", "ccbbaa23", "ccbbaa23", "ccbbaa23", "ccbbaa24", 
    "ccbbaa24", "ccbbaa24", "ccbbaa29", "ccbbaa29", "ccbbaa29", 
    "ccbbaa30", "ccbbaa30", "ccbbaa30", "ccbbaa31", "ccbbaa31", 
    "ccbbaa31", "ccbbaa33", "ccbbaa33", "ccbbaa33", "ccbbaa34", 
    "ccbbaa34", "ccbbaa34", "ccbbaa36", "ccbbaa36", "ccbbaa36", 
    "ccbbaa37", "ccbbaa37", "ccbbaa37", "ccbbaa38", "ccbbaa38", 
    "ccbbaa38", "ccbbaa41", "ccbbaa41", "ccbbaa41", "ccbbaa42", 
    "ccbbaa42", "ccbbaa42", "ccbbaa43", "ccbbaa43", "ccbbaa43", 
    "ccbbaa44", "ccbbaa44", "ccbbaa44", "ccbbaa45", "ccbbaa45", 
    "ccbbaa45", "ccbbaa47", "ccbbaa47", "ccbbaa47", "ccbbaa48", 
    "ccbbaa48", "ccbbaa48", "ccbbaa49", "ccbbaa49", "ccbbaa49", 
    "ccbbaa50", "ccbbaa50", "ccbbaa50", "ccbbaa51", "ccbbaa51", 
    "ccbbaa51", "ccbbaa52", "ccbbaa52", "ccbbaa52", "ccbbaa54", 
    "ccbbaa54", "ccbbaa54", "ccbbaa55", "ccbbaa55", "ccbbaa55", 
    "ccbbaa56", "ccbbaa56", "ccbbaa56", "ccbbaa57", "ccbbaa57", 
    "ccbbaa57", "ccbbaa58", "ccbbaa58", "ccbbaa58", "ccbbaa59", 
    "ccbbaa59", "ccbbaa59", "ccbbaa61", "ccbbaa61", "ccbbaa61", 
    "ccbbaa62", "ccbbaa62", "ccbbaa62", "ccbbaa63", "ccbbaa63", 
    "ccbbaa63", "ccbbaa64", "ccbbaa64", "ccbbaa64", "ccbbaa66", 
    "ccbbaa66", "ccbbaa66", "ccbbaa68", "ccbbaa68", "ccbbaa68", 
    "ccbbaa69", "ccbbaa69", "ccbbaa69", "ccbbaa70", "ccbbaa70", 
    "ccbbaa70", "ccbbaa71", "ccbbaa71", "ccbbaa71", "ccbbaa72", 
    "ccbbaa72", "ccbbaa72", "ccbbaa73", "ccbbaa73", "ccbbaa73", 
    "ccbbaa74", "ccbbaa74", "ccbbaa74", "ccbbaa75", "ccbbaa75", 
    "ccbbaa75", "ccbbaa76", "ccbbaa76", "ccbbaa76", "ccbbaa77", 
    "ccbbaa77", "ccbbaa77", "ccddaabb", "ccddaabb", "ccddaabb", 
    "ccddaacc", "ccddaacc", "ccddaacc", "ccddaa05", "ccddaa05", 
    "ccddaa05", "ccddaa09", "ccddaa09", "ccddaa09", "ccddaa14", 
    "ccddaa14", "ccddaa14", "ccddaa15", "ccddaa15", "ccddaa15", 
    "ccddaa17", "ccddaa17", "ccddaa17", "ccddaa20", "ccddaa20", 
    "ccddaa20", "ccddaa21", "ccddaa21", "ccddaa21", "ccddaa22", 
    "ccddaa22", "ccddaa22", "ccddaa26", "ccddaa26", "ccddaa26", 
    "ccddaa27", "ccddaa27", "ccddaa27", "ccddaa28", "ccddaa28", 
    "ccddaa28", "ccddaa29", "ccddaa29", "ccddaa29", "ccddaa30", 
    "ccddaa30", "ccddaa30", "ccddaa31", "ccddaa31", "ccddaa31", 
    "ccddaa32", "ccddaa32", "ccddaa32", "ccddaa33", "ccddaa33", 
    "ccddaa33", "ccddaa34", "ccddaa34", "ccddaa34", "ccddaa35", 
    "ccddaa35", "ccddaa35", "ccddaa36", "ccddaa36", "ccddaa36", 
    "ccddaa37", "ccddaa37", "ccddaa37", "ccddaa38", "ccddaa38", 
    "ccddaa38", "ccddaa39", "ccddaa39", "ccddaa39", "ccddaa40", 
    "ccddaa40", "ccddaa40", "ccddaa41", "ccddaa41", "ccddaa41", 
    "ccddaa42", "ccddaa42", "ccddaa42", "ccddaa43", "ccddaa43", 
    "ccddaa43", "ccddaa45", "ccddaa45", "ccddaa45", "ccddaa47", 
    "ccddaa47", "ccddaa47", "ccddaa48", "ccddaa48", "ccddaa48", 
    "ccddaa49", "ccddaa49", "ccddaa49", "ccddaa50", "ccddaa50", 
    "ccddaa50", "ccddaa51", "ccddaa51", "ccddaa51", "ccddaa52", 
    "ccddaa52", "ccddaa52", "ccccaabb", "ccccaabb", "ccccaabb", 
    "ccccaa04", "ccccaa04", "ccccaa04", "ccccaa05", "ccccaa05", 
    "ccccaa05", "ccccaa06", "ccccaa06", "ccccaa06", "ccccaa09", 
    "ccccaa09", "ccccaa09", "ccccaa12", "ccccaa12", "ccccaa12", 
    "ccccaa13", "ccccaa13", "ccccaa13", "ccccaa14", "ccccaa14", 
    "ccccaa14", "ccccaa15", "ccccaa15", "ccccaa15", "ccccaa16", 
    "ccccaa16", "ccccaa16", "ccccaa17", "ccccaa17", "ccccaa17", 
    "ccccaa18", "ccccaa18", "ccccaa18", "ccccaa19", "ccccaa19", 
    "ccccaa19", "ccccaa20", "ccccaa20", "ccccaa20", "ccccaa22", 
    "ccccaa22", "ccccaa22", "ccccaa23", "ccccaa23", "ccccaa23", 
    "ccccaa24", "ccccaa24", "ccccaa24", "ccccaa25", "ccccaa25", 
    "ccccaa25", "ccccaa26", "ccccaa26", "ccccaa26", "ccccaa27", 
    "ccccaa27", "ccccaa27", "ccccaa28", "ccccaa28", "ccccaa28", 
    "ccccaa30", "ccccaa30", "ccccaa30", "ccccaa31", "ccccaa31", 
    "ccccaa31", "ccccaa32", "ccccaa32", "ccccaa32"), Arm = structure(c(2L, 
    2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 
    2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
    2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 
    1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 
    2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
    2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
    1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L), levels = c("A", "B", "C"), class = "factor"), Visit = structure(c(1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L), levels = c("Month 6", "Month 12", "Month 20"), class = c("ordered", 
    "factor")), ODIPain = c(2L, 1L, 1L, 0L, 0L, 2L, 1L, 3L, 5L, 
    2L, 0L, 0L, 3L, 4L, 0L, 2L, 1L, 2L, 3L, 3L, 3L, 1L, 2L, 1L, 
    2L, 2L, 3L, 4L, 2L, 2L, 1L, 0L, 0L, 4L, 3L, 3L, 4L, 2L, 1L, 
    5L, 2L, 4L, 1L, 2L, 2L, 2L, 2L, 3L, 0L, 3L, 3L, 1L, 1L, 1L, 
    2L, 2L, 2L, 2L, 2L, 1L, 0L, 0L, 0L, 4L, 0L, 3L, 2L, 2L, 3L, 
    2L, 4L, 5L, 1L, 4L, 3L, 4L, 4L, 4L, 4L, 1L, 2L, 1L, 2L, 3L, 
    4L, 2L, 2L, 2L, 2L, 2L, 4L, 4L, 1L, 0L, 3L, 3L, 2L, 2L, 3L, 
    2L, 4L, 0L, 2L, 2L, 2L, 0L, 0L, 1L, 4L, 3L, 4L, 2L, 3L, 3L, 
    4L, 3L, 3L, 0L, 1L, 2L, 0L, 1L, 0L, 1L, 0L, 4L, 0L, 0L, 0L, 
    1L, 0L, 0L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 0L, 3L, 0L, 2L, 0L, 
    1L, 1L, 0L, 0L, 0L, 0L, 1L, 2L, 0L, 3L, 1L, 0L, 5L, 5L, 5L, 
    3L, 3L, 2L, 2L, 4L, 2L, 2L, 1L, 0L, 2L, 0L, 0L, 2L, 0L, 0L, 
    0L, 3L, 1L, 3L, 4L, 5L, 4L, 3L, 2L, 3L, 4L, 4L, 2L, 2L, 0L, 
    2L, 3L, 2L, 0L, 0L, 2L, 2L, 1L, 1L, 0L, 0L, 1L, 3L, 5L, 3L, 
    2L, 2L, 1L, 2L, 2L, 3L, 1L, 2L, 2L, 2L, 3L, 3L, 2L, 2L, 4L, 
    2L, 2L, 3L, 2L, 0L, 2L, 2L, 5L, 3L, 1L, 0L, 0L, 1L, 1L, 1L, 
    0L, 1L, 3L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 2L, 0L, 2L, 2L, 2L, 
    3L, 4L, 2L, 3L, 2L, 2L, 2L, 2L, 0L, 1L, 0L, 2L, 1L, 0L, 1L, 
    1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 2L, 0L, 1L, 3L, 0L, 0L, 0L, 
    2L, 1L, 1L, 0L, 0L, 0L, 2L, 2L, 1L, 0L, 0L, 0L, 2L, 3L, 2L, 
    3L, 5L, 5L, 3L, 3L, 2L, 2L, 4L, 2L, 2L, 2L, 2L, 0L, 0L, 0L, 
    1L, 1L, 2L, 2L, 1L, 2L, 2L, 4L, 3L, 2L, 4L, 2L, 2L, 2L, 2L, 
    4L, 4L, 3L, 1L, 2L, 1L, 1L, 0L, 1L)), row.names = c(56L, 
167L, 278L, 57L, 168L, 279L, 1L, 112L, 223L, 58L, 169L, 280L, 
2L, 113L, 224L, 59L, 170L, 281L, 3L, 114L, 225L, 60L, 171L, 282L, 
61L, 172L, 283L, 62L, 173L, 284L, 4L, 115L, 226L, 5L, 116L, 227L, 
6L, 117L, 228L, 63L, 174L, 285L, 64L, 175L, 286L, 7L, 118L, 229L, 
8L, 119L, 230L, 65L, 176L, 287L, 9L, 120L, 231L, 66L, 177L, 288L, 
67L, 178L, 289L, 10L, 121L, 232L, 68L, 179L, 290L, 11L, 122L, 
233L, 12L, 123L, 234L, 69L, 180L, 291L, 13L, 124L, 235L, 14L, 
125L, 236L, 70L, 181L, 292L, 71L, 182L, 293L, 72L, 183L, 294L, 
15L, 126L, 237L, 16L, 127L, 238L, 17L, 128L, 239L, 73L, 184L, 
295L, 18L, 129L, 240L, 19L, 130L, 241L, 74L, 185L, 296L, 20L, 
131L, 242L, 75L, 186L, 297L, 76L, 187L, 298L, 77L, 188L, 299L, 
21L, 132L, 243L, 22L, 133L, 244L, 78L, 189L, 300L, 23L, 134L, 
245L, 24L, 135L, 246L, 79L, 190L, 301L, 25L, 136L, 247L, 80L, 
191L, 302L, 26L, 137L, 248L, 81L, 192L, 303L, 27L, 138L, 249L, 
82L, 193L, 304L, 28L, 139L, 250L, 83L, 194L, 305L, 29L, 140L, 
251L, 84L, 195L, 306L, 85L, 196L, 307L, 86L, 197L, 308L, 87L, 
198L, 309L, 30L, 141L, 252L, 88L, 199L, 310L, 89L, 200L, 311L, 
31L, 142L, 253L, 32L, 143L, 254L, 33L, 144L, 255L, 90L, 201L, 
312L, 91L, 202L, 313L, 92L, 203L, 314L, 34L, 145L, 256L, 35L, 
146L, 257L, 36L, 147L, 258L, 37L, 148L, 259L, 38L, 149L, 260L, 
93L, 204L, 315L, 94L, 205L, 316L, 95L, 206L, 317L, 39L, 150L, 
261L, 96L, 207L, 318L, 97L, 208L, 319L, 40L, 151L, 262L, 98L, 
209L, 320L, 41L, 152L, 263L, 42L, 153L, 264L, 99L, 210L, 321L, 
43L, 154L, 265L, 100L, 211L, 322L, 44L, 155L, 266L, 45L, 156L, 
267L, 101L, 212L, 323L, 46L, 157L, 268L, 102L, 213L, 324L, 47L, 
158L, 269L, 48L, 159L, 270L, 103L, 214L, 325L, 49L, 160L, 271L, 
104L, 215L, 326L, 50L, 161L, 272L, 105L, 216L, 327L, 51L, 162L, 
273L, 52L, 163L, 274L, 106L, 217L, 328L, 107L, 218L, 329L, 53L, 
164L, 275L, 54L, 165L, 276L, 108L, 219L, 330L, 55L, 166L, 277L, 
109L, 220L, 331L, 110L, 221L, 332L, 111L, 222L, 333L), class = "data.frame"))
Generalized commented 2 months ago

OK, it was trivial...

(imp_mod_em_direct <- emmeans(as.mira(imp_mod), 
                       specs = ~ Visit | Arm, adjust="none",
                       vcov = pooled_coefficients$pooled_VC))  # <--- I can pass the pooled vcov matrix this way directly, when everything else is already established.