search

rvlenth / emmeans

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

emmeans seems to be losing the information about trend in multiple imputation setting #486

Closed Generalized closed 2 months ago

Generalized commented 2 months ago

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

> g <- glmgee(ODIPain ~ Visit * Arm,
+        id = PatientId,
+        data = d,
+        family = gaussian(link = "identity"),
+        corstr = "unstructured")
> coef(g)
             Estimates
(Intercept)     1.9030
Visit.L         0.2443
Visit.Q        -0.0965
ArmB           -0.2304
Visit.L:ArmB   -0.4589
Visit.Q:ArmB    0.0455

The model recognizes the ordered factor.

Now qdrg() (I need it for some reasons)

(qgrid <- qdrg(formula = ODIPain ~ Visit * Arm,
     data=d,
     coef = g$coefficients, 
     vcov = vcov(g, type = "bias-corrected")))

'emmGrid' object with variables:
    Visit = Month 6, Month 12, Month 20
    Arm = A, B
    rep.meas = multivariate response levels: 

> (qem <- emmeans(qgrid, specs = ~Visit | Arm))
Arm = A:
 Visit    emmean    SE  df asymp.LCL asymp.UCL
 Month 6    1.90 0.161 Inf      1.59      2.22
 Month 12   2.15 0.235 Inf      1.69      2.61
 Month 20   1.81 0.197 Inf      1.42      2.19

Arm = B:
 Visit    emmean    SE  df asymp.LCL asymp.UCL
 Month 6    1.67 0.137 Inf      1.40      1.94
 Month 12   1.46 0.192 Inf      1.08      1.83
 Month 20   1.62 0.157 Inf      1.31      1.93

Confidence level used: 0.95

OK, it works here:

> update(contrast(qem, "poly", max.degree=2, 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.096 0.133 Inf    -0.357     0.164  -0.727  0.4670
 quadratic   -0.585 0.316 Inf    -1.205     0.035  -1.850  0.0640

Arm = B:
 contrast  estimate    SE  df asymp.LCL asymp.UCL z.ratio p.value
 linear      -0.051 0.107 Inf    -0.260     0.158  -0.479  0.6320
 quadratic    0.378 0.308 Inf    -0.225     0.982   1.229  0.2190

Confidence level used: 0.95

Now trying obtain the same in multiple imputation manner. I cannot attach the data, can only describe the problem.

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

gives a list of analyses. Coefficients:

> geemmrm_imp_list_trend[[1]]$coefficients

(Intercept)              1.9030
Visit.L                  0.2443
Visit.Q                 -0.0965
ArmHD         -0.2304
Visit.L:ArmHD -0.4589
Visit.Q:ArmHD  0.0455

and

> vcov(geemmrm_imp_list_trend[[1]], type = "bias-corrected")
                        (Intercept)   Visit.L   Visit.Q ArmHD Visit.L:ArmHD Visit.Q:ArmHD 
(Intercept)                 0.02596  0.004650 -0.002323        -0.02596               -0.004650                0.002323
Visit.L                     0.00465  0.019965 -0.000644        -0.00465               -0.019965                0.000644
Visit.Q                    -0.00232 -0.000644  0.017600         0.00232                0.000644               -0.017600
ArmHD -0.02596 -0.004650  0.002323         0.04480                0.005531               -0.005062
Visit.L:ArmHD -0.00465 -0.019965  0.000644         0.00553                0.036137               -0.005332
Visit.Q:ArmHD 0.00232  0.000644 -0.017600        -0.00506               -0.005332                0.028957

Pooling it:

mmrm_gee_imp_grid_trend <- with(pool_estimates_for_qdrg(k=20, 
                                                  coefficients = lapply(geemmrm_imp_list_trend, function(analysis) analysis$coefficients),
                                                  varcov       = lapply(geemmrm_imp_list_trend, function(analysis) vcov(analysis, type = "bias-corrected"))),
                          qdrg(formula = ODIPain ~ Visit * Arm,
                               data=imp_data_long[[1]], # fake data in any case 
                               coef = pooled_coef, 
                               vcov = pooled_VC))

> mmrm_gee_imp_grid_trend
'emmGrid' object with variables:
    Visit = Month 6, Month 12, Month 20
    Arm = SoC, HD 

> mmrm_gee_imp_grid_trend@levels
$Visit
[1] "Month 6"  "Month 12" "Month 20"

$Arm
[1] "SoC"          "HD"

> mmrm_gee_imp_grid_trend@grid
     Visit          Arm .wgt.
1  Month 6          SoC    55
2 Month 12          SoC    55
3 Month 20          SoC    55
4  Month 6 HD 56
5 Month 12 HD 56
6 Month 20 HD 56

> mmrm_gee_imp_grid_trend@bhat
                  (Intercept)                 VisitMonth 12                 VisitMonth 20               ArmHD 
                        1.905                            NA                            NA                        -0.216 
VisitMonth 12:ArmHD VisitMonth 20:ArmHD 
                           NA                            NA

Maybe the bhat is the issue?

Now all estimates are the same:

> (mmrm_gee_imp_grid_trend_em <- emmeans(mmrm_gee_imp_grid_trend, 
+                                        specs = ~ Visit | Arm, adjust="none"))
Arm = SoC:
 Visit    emmean    SE  df asymp.LCL asymp.UCL
 Month 6    1.90 0.163 Inf      1.59      2.22
 Month 12   1.90 0.163 Inf      1.59      2.22
 Month 20   1.90 0.163 Inf      1.59      2.22

Arm = HD :
 Visit    emmean    SE  df asymp.LCL asymp.UCL
 Month 6    1.69 0.142 Inf      1.41      1.97
 Month 12   1.69 0.142 Inf      1.41      1.97
 Month 20   1.69 0.142 Inf      1.41      1.97

Confidence level used: 0.95

and finally

> update(contrast(mmrm_gee_imp_grid_trend_em, "poly", by = "Arm"), 
+        adjust="none", level = 0.95, infer = c(TRUE, TRUE))
Arm = SoC:
 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 = HD :
 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

Do you have any idea, form this limited data, what could cause that?

I can provide more details about the grid object if you request, but providing all the data will be challenging...

Data for the simple, working case.:

d <- > dput(d)
structure(list(PatientId = c("xxayyy04", "xxayyy04", "xxayyy04", 
"xxayyy07", "xxayyy07", "xxayyy07", "xxayyy08", "xxayyy08", "xxayyy08", 
"xxayyy10", "xxayyy10", "xxayyy10", "xxayyy12", "xxayyy12", "xxayyy12", 
"xxayyy14", "xxayyy14", "xxayyy14", "xxayyy17", "xxayyy17", "xxayyy17", 
"xxayyy18", "xxayyy18", "xxayyy18", "xxayyy19", "xxayyy19", "xxayyy19", 
"xxayyy21", "xxayyy21", "xxayyy21", "xxayyy23", "xxayyy23", "xxayyy23", 
"xxayyy24", "xxayyy24", "xxayyy24", "xxayyy29", "xxayyy29", "xxayyy29", 
"xxay0xx0", "xxay0xx0", "xxay0xx0", "xxay0xx1", "xxay0xx1", "xxay0xx1", 
"xxay0xx3", "xxay0xx3", "xxay0xx3", "xxay0xx4", "xxay0xx4", "xxay0xx4", 
"xxay0xx6", "xxay0xx6", "xxay0xx6", "xxay0xx7", "xxay0xx7", "xxay0xx7", 
"xxay0xx8", "xxay0xx8", "xxay0xx8", "xxayyy41", "xxayyy41", "xxayyy41", 
"xxayyy42", "xxayyy42", "xxayyy42", "xxayyy43", "xxayyy43", "xxayyy43", 
"xxayyy44", "xxayyy44", "xxayyy44", "xxayyy45", "xxayyy45", "xxayyy45", 
"xxayyy47", "xxayyy47", "xxayyy47", "xxayyy48", "xxayyy48", "xxayyy48", 
"xxayyy49", "xxayyy49", "xxayyy49", "xxayyy50", "xxayyy50", "xxayyy50", 
"xxayyy51", "xxayyy51", "xxayyy51", "xxayyy52", "xxayyy52", "xxayyy52", 
"xxayyy54", "xxayyy54", "xxayyy54", "xxayyy55", "xxayyy55", "xxayyy55", 
"xxayyy56", "xxayyy56", "xxayyy56", "xxayyy57", "xxayyy57", "xxayyy57", 
"xxayyy58", "xxayyy58", "xxayyy58", "xxayyy59", "xxayyy59", "xxayyy59", 
"xxayyy61", "xxayyy61", "xxayyy61", "xxayyy62", "xxayyy62", "xxayyy62", 
"xxayyy63", "xxayyy63", "xxayyy63", "xxayyy64", "xxayyy64", "xxayyy64", 
"xxayyy66", "xxayyy66", "xxayyy66", "xxayyy68", "xxayyy68", "xxayyy68", 
"xxayyy69", "xxayyy69", "xxayyy69", "xxayyy70", "xxayyy70", "xxayyy70", 
"xxayyy71", "xxayyy71", "xxayyy71", "xxayyy72", "xxayyy72", "xxayyy72", 
"xxayyy73", "xxayyy73", "xxayyy73", "xxayyy74", "xxayyy74", "xxayyy74", 
"xxayyy75", "xxayyy75", "xxayyy75", "xxayyy76", "xxayyy76", "xxayyy76", 
"xxayyy77", "xxayyy77", "xxayyy77", "xxbyyyay", "xxbyyyay", "xxbyyyay", 
"xxbyyyxx", "xxbyyyxx", "xxbyyyxx", "xxbyyy05", "xxbyyy05", "xxbyyy05", 
"xxbyyy09", "xxbyyy09", "xxbyyy09", "xxbyyy14", "xxbyyy14", "xxbyyy14", 
"xxbyyy15", "xxbyyy15", "xxbyyy15", "xxbyyy17", "xxbyyy17", "xxbyyy17", 
"xxbyyy20", "xxbyyy20", "xxbyyy20", "xxbyyy21", "xxbyyy21", "xxbyyy21", 
"xxbyyy22", "xxbyyy22", "xxbyyy22", "xxbyyy26", "xxbyyy26", "xxbyyy26", 
"xxbyyy27", "xxbyyy27", "xxbyyy27", "xxbyyy28", "xxbyyy28", "xxbyyy28", 
"xxbyyy29", "xxbyyy29", "xxbyyy29", "xxby0xx0", "xxby0xx0", "xxby0xx0", 
"xxby0xx1", "xxby0xx1", "xxby0xx1", "xxby0xx2", "xxby0xx2", "xxby0xx2", 
"xxby0xx3", "xxby0xx3", "xxby0xx3", "xxby0xx4", "xxby0xx4", "xxby0xx4", 
"xxby0xx5", "xxby0xx5", "xxby0xx5", "xxby0xx6", "xxby0xx6", "xxby0xx6", 
"xxby0xx7", "xxby0xx7", "xxby0xx7", "xxby0xx8", "xxby0xx8", "xxby0xx8", 
"xxby0xx9", "xxby0xx9", "xxby0xx9", "xxbyyy40", "xxbyyy40", "xxbyyy40", 
"xxbyyy41", "xxbyyy41", "xxbyyy41", "xxbyyy42", "xxbyyy42", "xxbyyy42", 
"xxbyyy43", "xxbyyy43", "xxbyyy43", "xxbyyy45", "xxbyyy45", "xxbyyy45", 
"xxbyyy47", "xxbyyy47", "xxbyyy47", "xxbyyy48", "xxbyyy48", "xxbyyy48", 
"xxbyyy49", "xxbyyy49", "xxbyyy49", "xxbyyy50", "xxbyyy50", "xxbyyy50", 
"xxbyyy51", "xxbyyy51", "xxbyyy51", "xxbyyy52", "xxbyyy52", "xxbyyy52", 
"xxxxyyay", "xxxxyyay", "xxxxyyay", "xxxxyy04", "xxxxyy04", "xxxxyy04", 
"xxxxyy05", "xxxxyy05", "xxxxyy05", "xxxxyy06", "xxxxyy06", "xxxxyy06", 
"xxxxyy09", "xxxxyy09", "xxxxyy09", "xxxxyy12", "xxxxyy12", "xxxxyy12", 
"xxxxyy13", "xxxxyy13", "xxxxyy13", "xxxxyy14", "xxxxyy14", "xxxxyy14", 
"xxxxyy15", "xxxxyy15", "xxxxyy15", "xxxxyy16", "xxxxyy16", "xxxxyy16", 
"xxxxyy17", "xxxxyy17", "xxxxyy17", "xxxxyy18", "xxxxyy18", "xxxxyy18", 
"xxxxyy19", "xxxxyy19", "xxxxyy19", "xxxxyy20", "xxxxyy20", "xxxxyy20", 
"xxxxyy22", "xxxxyy22", "xxxxyy22", "xxxxyy23", "xxxxyy23", "xxxxyy23", 
"xxxxyy24", "xxxxyy24", "xxxxyy24", "xxxxyy25", "xxxxyy25", "xxxxyy25", 
"xxxxyy26", "xxxxyy26", "xxxxyy26", "xxxxyy27", "xxxxyy27", "xxxxyy27", 
"xxxxyy28", "xxxxyy28", "xxxxyy28", "xxxx0xx0", "xxxx0xx0", "xxxx0xx0", 
"xxxx0xx1", "xxxx0xx1", "xxxx0xx1", "xxxx0xx2", "xxxx0xx2", "xxxx0xx2"
), 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), 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")), 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")), 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

Will try to provide limited imputed dataset later.