1
$\begingroup$

Background

I am writing a project on a big multicenter RCT study, where subjects are following different dietary patterns for 2 years. I have access to their dietary intakes and I am grouping them in quartiles of protein intake.

My aim is to test whether there are any differences in effects on outcomes (continuous variables like body weight, blood parameters, etc.) between these groups.

My Approach

My approach has been to fit the full model with lme() e.g.

modWeight <- lme(wght~protein_quan*factor(week)+w0+energy_kj+age+gen+eth+smoking_status+medication+baecke_mean+alcohol,
                    random=~1|site/subject_id,method="ML",
                    corr = corExp(form = ~ week | site/subject_id, nugget=T),
                    data_wght, na.action=na.exclude)

Where: - protein_quan is a factor of the protein intake quartile groups - factor(week) is a factor of the weeks, w0 is a continuous covariate of baseline weight - energy_kj is a continuous covariate of energy intake - age is a continuous covariate of subject age - gen is a factor of the subject sex group - eth is a factor of the subject ethnic group - smoking_status is a factor of the smoking status - medication is a factor of medication intake - baecke_mean is a continuos covariate of Physical Activity

I then proceed to reduce the fitted model to the most until I reach the one that is parsimonious following these steps:

  1. Test if the model can be reduced to the continuous variable of week (time) with anova(modContinuous,modFactor)i
  2. Analyse distribution of model, transform response if necessary (log(),sqrt(), or boxcox)
  3. Test to see if I can reduce for serial correlation (none vs corGaus / corExp)
  4. Reduce the random effects
  5. Reduce the fixed effects
  6. Final model:

``

modfitWeight <-lme(wght_cox~protein_quan+factor(week)+w0+fd_energy_kj+smoking_status,
                    random=~1|subject_id,method="REML",
                    corr = corGaus(form = ~ week | subject_id, nugget=T),
                    data=data_wght_trans,na.action=na.exclude)

Where the response variable has been transformed with boxcox.

Find my data here.

I plot the mean fitted response variables grouped by protein group and week, with ggplot:

data_wght_cmplt_trans_kgbwxenergy %>%
  mutate(pred_wght = (lambda_wght*fitted(modfitWeight)+1)^(1/lambda_wght)) %>%
  group_by(week,protein_quan)%>%
  mutate(pred_w0=pred_wght[week==min(week)],
         pred_wght_mean=mean(pred_wght),
         pred_wght_mean_upper = mean(pred_wght) + qnorm(0.975) * (sd(pred_wght) / sqrt(sum(!is.na(pred_wght)))),
         pred_wght_mean_lower = mean(pred_wght) - qnorm(0.975) * (sd(pred_wght) / sqrt(sum(!is.na(pred_wght)))),
         pred_wght_mean_w0=mean(pred_wght-pred_w0)) %>%
  ungroup() %>%
  ggplot(aes(x = week, y=pred_wght_mean,color=protein_quan)) + 
  geom_line(size = 1) +
  geom_errorbar(aes(ymin = pred_wght_mean_lower, ymax = pred_wght_mean_upper)) + 
  theme_pubr() 

ggplot of adjusted weight means

To analyse if there are any significant differences between the protein intake groups or between different weeks, I first take a look at the summary of the model:

summary(modfitWeight)

Linear mixed-effects model fit by REML
 Data: data_wght_cmplt_trans_kgbwxenergy 
       AIC      BIC   logLik
  191.4532 270.5222 -74.7266

Random effects:
 Formula: ~1 | subject_id
         (Intercept)  Residual
StdDev: 0.0002939717 0.7757884

Correlation Structure: Gaussian spatial correlation
 Formula: ~week | subject_id 
 Parameter estimate(s):
      range      nugget 
44.71884151  0.03657875 
Fixed effects: wght_cox ~ protein_quan + factor(week) + w0 + smoking_status 
                              Value Std.Error  DF   t-value p-value
(Intercept)                9.869963 0.8066643 297 12.235527  0.0000
protein_quan2             -0.508055 0.2880310  22 -1.763891  0.0916
protein_quan3             -0.707052 0.3104516  22 -2.277495  0.0328
protein_quan4             -0.690427 0.3041147  22 -2.270286  0.0333
factor(week)10            -0.044499 0.0406847 297 -1.093748  0.2750
factor(week)12            -0.122841 0.0436181 297 -2.816294  0.0052
factor(week)16            -0.214534 0.0536375 297 -3.999701  0.0001
factor(week)20            -0.259108 0.0667081 297 -3.884200  0.0001
factor(week)26            -0.431458 0.0881326 297 -4.895552  0.0000
factor(week)32            -0.322749 0.1092590 297 -2.953982  0.0034
factor(week)44            -0.260876 0.1460344 297 -1.786400  0.0751
factor(week)52            -0.341910 0.1651032 297 -2.070888  0.0392
factor(week)64            -0.149894 0.1853560 297 -0.808681  0.4193
factor(week)78            -0.064264 0.1985349 297 -0.323689  0.7464
factor(week)104            0.125741 0.2063403 297  0.609386  0.5427
w0                         0.111063 0.0082458  22 13.469155  0.0000
smoking_statusDaily       -1.434396 0.5635267  22 -2.545392  0.0184

The average subject weight for week 8 (baseline) and protein group 1 is 96.4:

data_wght %>%
    mutate(pred_wght = (lambda_wght*fitted(modfitWeight)+1)^(1/lambda_wght)) %>%
    filter(week==8,fd_protein_quan_kgbw_adj=="1") %>%
    mutate(mean_weight_group=mean(pred_wght)) %>%
    dplyr::select(week, fd_protein_quan_kgbw_adj, mean_weight_group)

# A tibble: 7 x 3
   week protein_quan      mean_weight_group
  <dbl> <fct>             <dbl>
1     8 1                 96.4

By back transforming the boxcox response I get an adjusted response for week 8 and protein group 1 of:

> (lambda_cmplt_wght_kgbwxenergy*(9.869963+96.4*0.111063)+1)^(1/lambda_cmplt_wght_kgbwxenergy)
[1] 98.44476

Now I want to compare contrasts and see if there are any significant differences in my modelled groups. I try with emmeans and get a summary of my back-transformed boxcox means:

emmeans.wght <- emmeans(modfitWeight, ~ protein_quan|week )
emmeans.wght_cox <- update(emmeans.wght,tran=make.tran("boxcox", lambda_wght))
summary(emmeans.wght_cox, type="response")

 protein_quan = 1:
 week    response       SE df lower.CL upper.CL
 8       82.83042 2.740448 22 77.23589 88.60220
 10      82.49940 2.735465 22 76.91523 88.26088
 12      81.91808 2.726686 22 76.35218 87.66141
 16      81.24007 2.716405 22 75.69556 86.96215
 20      80.91140 2.711404 22 75.37730 86.62314
 26      79.64627 2.692051 22 74.15243 85.31800
 32      80.44319 2.704261 22 74.92395 86.14016
 44      80.89838 2.711206 22 75.36469 86.60971
 52      80.30246 2.702110 22 74.78770 85.99499
 64      81.71778 2.723654 22 76.15819 87.45484
 78      82.35256 2.733251 22 76.77301 88.10947
 104     83.76905 2.754521 22 78.14524 89.56993

 protein_quan = 2:
 week    response       SE df lower.CL upper.CL
 8       79.08692 2.629174 22 73.71996 84.62475
 10      78.76279 2.624271 22 73.40603 84.29048
 12      78.19361 2.615633 22 72.85483 83.70345
 16      77.52983 2.605517 22 72.21209 83.01875
 20      77.20808 2.600597 22 71.90058 82.68683
 26      75.96972 2.581553 22 70.70184 81.40910
 32      76.74975 2.593568 22 71.45687 82.21397
 44      77.19533 2.600401 22 71.88823 82.67368
 52      76.61200 2.591451 22 71.32352 82.07184
 64      77.99751 2.612650 22 72.66493 83.50118
 78      78.61902 2.622092 22 73.26679 84.14221
 104     80.00608 2.643020 22 74.61032 85.57254

 protein_quan = 3:
 week    response       SE df lower.CL upper.CL
 8       77.64213 2.198794 22 73.14308 82.26291
 10      77.32071 2.194652 22 72.83028 81.93292
 12      76.75630 2.187355 22 72.28104 81.35342
 16      76.09810 2.178809 22 71.64062 80.67755
 20      75.77907 2.174653 22 71.33023 80.34992
 26      74.55123 2.158565 22 70.13584 79.08880
 32      75.32462 2.168715 22 70.88813 79.88319
 44      75.76643 2.174488 22 71.31793 80.33693
 52      75.18804 2.166927 22 70.75526 79.74291
 64      76.56184 2.184835 22 72.09183 81.15375
 78      77.17814 2.192811 22 72.69153 81.78654
 104     78.55365 2.210490 22 74.03028 83.19862

 protein_quan = 4:
 week    response       SE df lower.CL upper.CL
 8       77.76236 2.631738 22 72.39172 83.30713
 10      77.44072 2.626784 22 72.08038 82.97524
 12      76.87591 2.618058 22 71.53373 82.39240
 16      76.21725 2.607838 22 70.89633 81.71261
 20      75.89799 2.602867 22 70.58742 81.38307
 26      74.66927 2.583628 22 69.39873 80.11458
 32      75.44321 2.595767 22 70.14742 80.91362
 44      75.88534 2.602670 22 70.57518 81.37002
 52      75.30653 2.593628 22 70.01519 80.77252
 64      76.68132 2.615044 22 71.34541 82.19157
 78      77.29805 2.624583 22 71.94229 82.82803
 104     78.67452 2.645725 22 73.27478 84.24821

Results are averaged over the levels of: smoking_status 
Degrees-of-freedom method: containment 
Confidence level used: 0.95 
Intervals are back-transformed from the Box-Cox (lambda = 0.545) scale 

And by plotting the estimates with emmip:

emmip(emmeans.wght_cox, fd_protein_quan_kgbw_adj ~ week,CIs=TRUE)+theme_pubr()

emmip plot

I see that the emmeans estimates have the same trends, but are not in complete alignment in, respect to the estimate sizes, with the lme summary output, and the confidence interval bars are also different. This is when i cried my first tear

Then I tried the pairs function to see if there are any significant differences:

    protein_quan = 1:
 contrast     estimate         SE  df t.ratio p.value
 8 - 10    0.044498835 0.04068474 297   1.094  0.9948
 8 - 12    0.122841338 0.04361807 297   2.816  0.1786
 8 - 16    0.214534083 0.05363752 297   4.000  0.0045
 8 - 20    0.259107766 0.06670815 297   3.884  0.0070
 8 - 26    0.431457663 0.08813258 297   4.896  0.0001
 8 - 32    0.322749062 0.10925899 297   2.954  0.1279
 8 - 44    0.260875844 0.14603437 297   1.786  0.8243
 8 - 52    0.341910138 0.16510318 297   2.071  0.6439
 8 - 64    0.149893859 0.18535605 297   0.809  0.9997
 8 - 78    0.064263508 0.19853493 297   0.324  1.0000
 8 - 104  -0.125740959 0.20634035 297  -0.609  1.0000
 10 - 12   0.078342504 0.04068474 297   1.926  0.7423
 10 - 16   0.170035249 0.04807703 297   3.537  0.0234
 10 - 20   0.214608932 0.05993849 297   3.580  0.0202
 10 - 26   0.386958829 0.08092640 297   4.782  0.0002
 10 - 32   0.278250227 0.10235049 297   2.719  0.2224
 10 - 44   0.216377009 0.14055285 297   1.539  0.9283
 10 - 52   0.297411304 0.16076764 297   1.850  0.7887
 10 - 64   0.105395024 0.18262804 297   0.577  1.0000
 10 - 78   0.019764674 0.19719872 297   0.100  1.0000
 10 - 104 -0.170239794 0.20613091 297  -0.826  0.9996
 12 - 16   0.091692745 0.04361807 297   2.102  0.6216
 12 - 20   0.136266428 0.05363752 297   2.541  0.3192
 12 - 26   0.308616325 0.07375087 297   4.185  0.0022
 12 - 32   0.199907724 0.09529330 297   2.098  0.6248
 12 - 44   0.138034505 0.13479517 297   1.024  0.9971
 12 - 52   0.219068800 0.15614450 297   1.403  0.9623
 12 - 64   0.027052521 0.17965340 297   0.151  1.0000
 12 - 78  -0.058577830 0.19570155 297  -0.299  1.0000
 12 - 104 -0.248582297 0.20588322 297  -1.207  0.9881
 16 - 20   0.044573683 0.04361807 297   1.022  0.9971
 16 - 26   0.216923580 0.05993849 297   3.619  0.0177
 16 - 32   0.108214978 0.08092640 297   1.337  0.9735
 16 - 44   0.046341760 0.12249476 297   0.378  1.0000
 16 - 52   0.127376055 0.14603437 297   0.872  0.9993
 16 - 64  -0.064640225 0.17292481 297  -0.374  1.0000
 16 - 78  -0.150270575 0.19217335 297  -0.782  0.9998
 16 - 104 -0.340275043 0.20524946 297  -1.658  0.8855
 20 - 26   0.172349897 0.04807703 297   3.585  0.0199
 20 - 32   0.063641296 0.06670815 297   0.954  0.9984
 20 - 44   0.001768077 0.10925899 297   0.016  1.0000
 20 - 52   0.082802372 0.13479517 297   0.614  1.0000
 20 - 64  -0.109213908 0.16510318 297  -0.661  1.0000
 20 - 78  -0.194844258 0.18784840 297  -1.037  0.9967
 20 - 104 -0.384848725 0.20438557 297  -1.883  0.7690
 26 - 32  -0.108708602 0.04807703 297  -2.261  0.5068
 26 - 44  -0.170581820 0.08813258 297  -1.936  0.7359
 26 - 52  -0.089547525 0.11598318 297  -0.772  0.9998
 26 - 64  -0.281563805 0.15123292 297  -1.862  0.7817
 26 - 78  -0.367194155 0.17965340 297  -2.044  0.6629
 26 - 104 -0.557198623 0.20251567 297  -2.751  0.2069
 32 - 44  -0.061873218 0.06670815 297  -0.928  0.9988
 32 - 52   0.019161077 0.09529330 297   0.201  1.0000
 32 - 64  -0.172855203 0.13479517 297  -1.282  0.9808
 32 - 78  -0.258485554 0.16915399 297  -1.528  0.9317
 32 - 104 -0.448490021 0.19972270 297  -2.246  0.5180
 44 - 52   0.081034295 0.05363752 297   1.511  0.9367
 44 - 64  -0.110981985 0.09529330 297  -1.165  0.9911
 44 - 78  -0.196612335 0.14055285 297  -1.399  0.9631
 44 - 104 -0.386616803 0.19011677 297  -2.034  0.6701
 52 - 64  -0.192016280 0.06670815 297  -2.878  0.1542
 52 - 78  -0.277646630 0.11598318 297  -2.394  0.4135
 52 - 104 -0.467651098 0.17965340 297  -2.603  0.2828
 64 - 78  -0.085630351 0.07375087 297  -1.161  0.9914
 64 - 104 -0.275634818 0.15614450 297  -1.765  0.8354
 78 - 104 -0.190004467 0.11598318 297  -1.638  0.8935

fd_protein_quan_kgbw_adj = 2:
 contrast     estimate         SE  df t.ratio p.value
 8 - 10    0.044498835 0.04068474 297   1.094  0.9948
 8 - 12    0.122841338 0.04361807 297   2.816  0.1786
 8 - 16    0.214534083 0.05363752 297   4.000  0.0045
 8 - 20    0.259107766 0.06670815 297   3.884  0.0070
 8 - 26    0.431457663 0.08813258 297   4.896  0.0001
 8 - 32    0.322749062 0.10925899 297   2.954  0.1279
 8 - 44    0.260875844 0.14603437 297   1.786  0.8243
 8 - 52    0.341910138 0.16510318 297   2.071  0.6439
 8 - 64    0.149893859 0.18535605 297   0.809  0.9997
 8 - 78    0.064263508 0.19853493 297   0.324  1.0000
 8 - 104  -0.125740959 0.20634035 297  -0.609  1.0000
 10 - 12   0.078342504 0.04068474 297   1.926  0.7423
 10 - 16   0.170035249 0.04807703 297   3.537  0.0234
 10 - 20   0.214608932 0.05993849 297   3.580  0.0202
 10 - 26   0.386958829 0.08092640 297   4.782  0.0002
 10 - 32   0.278250227 0.10235049 297   2.719  0.2224
 10 - 44   0.216377009 0.14055285 297   1.539  0.9283
 10 - 52   0.297411304 0.16076764 297   1.850  0.7887
 10 - 64   0.105395024 0.18262804 297   0.577  1.0000
 10 - 78   0.019764674 0.19719872 297   0.100  1.0000
 10 - 104 -0.170239794 0.20613091 297  -0.826  0.9996
 12 - 16   0.091692745 0.04361807 297   2.102  0.6216
 12 - 20   0.136266428 0.05363752 297   2.541  0.3192
 12 - 26   0.308616325 0.07375087 297   4.185  0.0022
 12 - 32   0.199907724 0.09529330 297   2.098  0.6248
 12 - 44   0.138034505 0.13479517 297   1.024  0.9971
 12 - 52   0.219068800 0.15614450 297   1.403  0.9623
 12 - 64   0.027052521 0.17965340 297   0.151  1.0000
 12 - 78  -0.058577830 0.19570155 297  -0.299  1.0000
 12 - 104 -0.248582297 0.20588322 297  -1.207  0.9881
 16 - 20   0.044573683 0.04361807 297   1.022  0.9971
 16 - 26   0.216923580 0.05993849 297   3.619  0.0177
 16 - 32   0.108214978 0.08092640 297   1.337  0.9735
 16 - 44   0.046341760 0.12249476 297   0.378  1.0000
 16 - 52   0.127376055 0.14603437 297   0.872  0.9993
 16 - 64  -0.064640225 0.17292481 297  -0.374  1.0000
 16 - 78  -0.150270575 0.19217335 297  -0.782  0.9998
 16 - 104 -0.340275043 0.20524946 297  -1.658  0.8855
 20 - 26   0.172349897 0.04807703 297   3.585  0.0199
 20 - 32   0.063641296 0.06670815 297   0.954  0.9984
 20 - 44   0.001768077 0.10925899 297   0.016  1.0000
 20 - 52   0.082802372 0.13479517 297   0.614  1.0000
 20 - 64  -0.109213908 0.16510318 297  -0.661  1.0000
 20 - 78  -0.194844258 0.18784840 297  -1.037  0.9967
 20 - 104 -0.384848725 0.20438557 297  -1.883  0.7690
 26 - 32  -0.108708602 0.04807703 297  -2.261  0.5068
 26 - 44  -0.170581820 0.08813258 297  -1.936  0.7359
 26 - 52  -0.089547525 0.11598318 297  -0.772  0.9998
 26 - 64  -0.281563805 0.15123292 297  -1.862  0.7817
 26 - 78  -0.367194155 0.17965340 297  -2.044  0.6629
 26 - 104 -0.557198623 0.20251567 297  -2.751  0.2069
 32 - 44  -0.061873218 0.06670815 297  -0.928  0.9988
 32 - 52   0.019161077 0.09529330 297   0.201  1.0000
 32 - 64  -0.172855203 0.13479517 297  -1.282  0.9808
 32 - 78  -0.258485554 0.16915399 297  -1.528  0.9317
 32 - 104 -0.448490021 0.19972270 297  -2.246  0.5180
 44 - 52   0.081034295 0.05363752 297   1.511  0.9367
 44 - 64  -0.110981985 0.09529330 297  -1.165  0.9911
 44 - 78  -0.196612335 0.14055285 297  -1.399  0.9631
 44 - 104 -0.386616803 0.19011677 297  -2.034  0.6701
 52 - 64  -0.192016280 0.06670815 297  -2.878  0.1542
 52 - 78  -0.277646630 0.11598318 297  -2.394  0.4135
 52 - 104 -0.467651098 0.17965340 297  -2.603  0.2828
 64 - 78  -0.085630351 0.07375087 297  -1.161  0.9914
 64 - 104 -0.275634818 0.15614450 297  -1.765  0.8354
 78 - 104 -0.190004467 0.11598318 297  -1.638  0.8935

fd_protein_quan_kgbw_adj = 3:
 contrast     estimate         SE  df t.ratio p.value
 8 - 10    0.044498835 0.04068474 297   1.094  0.9948
 8 - 12    0.122841338 0.04361807 297   2.816  0.1786
 8 - 16    0.214534083 0.05363752 297   4.000  0.0045
 8 - 20    0.259107766 0.06670815 297   3.884  0.0070
 8 - 26    0.431457663 0.08813258 297   4.896  0.0001
 8 - 32    0.322749062 0.10925899 297   2.954  0.1279
 8 - 44    0.260875844 0.14603437 297   1.786  0.8243
 8 - 52    0.341910138 0.16510318 297   2.071  0.6439
 8 - 64    0.149893859 0.18535605 297   0.809  0.9997
 8 - 78    0.064263508 0.19853493 297   0.324  1.0000
 8 - 104  -0.125740959 0.20634035 297  -0.609  1.0000
 10 - 12   0.078342504 0.04068474 297   1.926  0.7423
 10 - 16   0.170035249 0.04807703 297   3.537  0.0234
 10 - 20   0.214608932 0.05993849 297   3.580  0.0202
 10 - 26   0.386958829 0.08092640 297   4.782  0.0002
 10 - 32   0.278250227 0.10235049 297   2.719  0.2224
 10 - 44   0.216377009 0.14055285 297   1.539  0.9283
 10 - 52   0.297411304 0.16076764 297   1.850  0.7887
 10 - 64   0.105395024 0.18262804 297   0.577  1.0000
 10 - 78   0.019764674 0.19719872 297   0.100  1.0000
 10 - 104 -0.170239794 0.20613091 297  -0.826  0.9996
 12 - 16   0.091692745 0.04361807 297   2.102  0.6216
 12 - 20   0.136266428 0.05363752 297   2.541  0.3192
 12 - 26   0.308616325 0.07375087 297   4.185  0.0022
 12 - 32   0.199907724 0.09529330 297   2.098  0.6248
 12 - 44   0.138034505 0.13479517 297   1.024  0.9971
 12 - 52   0.219068800 0.15614450 297   1.403  0.9623
 12 - 64   0.027052521 0.17965340 297   0.151  1.0000
 12 - 78  -0.058577830 0.19570155 297  -0.299  1.0000
 12 - 104 -0.248582297 0.20588322 297  -1.207  0.9881
 16 - 20   0.044573683 0.04361807 297   1.022  0.9971
 16 - 26   0.216923580 0.05993849 297   3.619  0.0177
 16 - 32   0.108214978 0.08092640 297   1.337  0.9735
 16 - 44   0.046341760 0.12249476 297   0.378  1.0000
 16 - 52   0.127376055 0.14603437 297   0.872  0.9993
 16 - 64  -0.064640225 0.17292481 297  -0.374  1.0000
 16 - 78  -0.150270575 0.19217335 297  -0.782  0.9998
 16 - 104 -0.340275043 0.20524946 297  -1.658  0.8855
 20 - 26   0.172349897 0.04807703 297   3.585  0.0199
 20 - 32   0.063641296 0.06670815 297   0.954  0.9984
 20 - 44   0.001768077 0.10925899 297   0.016  1.0000
 20 - 52   0.082802372 0.13479517 297   0.614  1.0000
 20 - 64  -0.109213908 0.16510318 297  -0.661  1.0000
 20 - 78  -0.194844258 0.18784840 297  -1.037  0.9967
 20 - 104 -0.384848725 0.20438557 297  -1.883  0.7690
 26 - 32  -0.108708602 0.04807703 297  -2.261  0.5068
 26 - 44  -0.170581820 0.08813258 297  -1.936  0.7359
 26 - 52  -0.089547525 0.11598318 297  -0.772  0.9998
 26 - 64  -0.281563805 0.15123292 297  -1.862  0.7817
 26 - 78  -0.367194155 0.17965340 297  -2.044  0.6629
 26 - 104 -0.557198623 0.20251567 297  -2.751  0.2069
 32 - 44  -0.061873218 0.06670815 297  -0.928  0.9988
 32 - 52   0.019161077 0.09529330 297   0.201  1.0000
 32 - 64  -0.172855203 0.13479517 297  -1.282  0.9808
 32 - 78  -0.258485554 0.16915399 297  -1.528  0.9317
 32 - 104 -0.448490021 0.19972270 297  -2.246  0.5180
 44 - 52   0.081034295 0.05363752 297   1.511  0.9367
 44 - 64  -0.110981985 0.09529330 297  -1.165  0.9911
 44 - 78  -0.196612335 0.14055285 297  -1.399  0.9631
 44 - 104 -0.386616803 0.19011677 297  -2.034  0.6701
 52 - 64  -0.192016280 0.06670815 297  -2.878  0.1542
 52 - 78  -0.277646630 0.11598318 297  -2.394  0.4135
 52 - 104 -0.467651098 0.17965340 297  -2.603  0.2828
 64 - 78  -0.085630351 0.07375087 297  -1.161  0.9914
 64 - 104 -0.275634818 0.15614450 297  -1.765  0.8354
 78 - 104 -0.190004467 0.11598318 297  -1.638  0.8935

fd_protein_quan_kgbw_adj = 4:
 contrast     estimate         SE  df t.ratio p.value
 8 - 10    0.044498835 0.04068474 297   1.094  0.9948
 8 - 12    0.122841338 0.04361807 297   2.816  0.1786
 8 - 16    0.214534083 0.05363752 297   4.000  0.0045
 8 - 20    0.259107766 0.06670815 297   3.884  0.0070
 8 - 26    0.431457663 0.08813258 297   4.896  0.0001
 8 - 32    0.322749062 0.10925899 297   2.954  0.1279
 8 - 44    0.260875844 0.14603437 297   1.786  0.8243
 8 - 52    0.341910138 0.16510318 297   2.071  0.6439
 8 - 64    0.149893859 0.18535605 297   0.809  0.9997
 8 - 78    0.064263508 0.19853493 297   0.324  1.0000
 8 - 104  -0.125740959 0.20634035 297  -0.609  1.0000
 10 - 12   0.078342504 0.04068474 297   1.926  0.7423
 10 - 16   0.170035249 0.04807703 297   3.537  0.0234
 10 - 20   0.214608932 0.05993849 297   3.580  0.0202
 10 - 26   0.386958829 0.08092640 297   4.782  0.0002
 10 - 32   0.278250227 0.10235049 297   2.719  0.2224
 10 - 44   0.216377009 0.14055285 297   1.539  0.9283
 10 - 52   0.297411304 0.16076764 297   1.850  0.7887
 10 - 64   0.105395024 0.18262804 297   0.577  1.0000
 10 - 78   0.019764674 0.19719872 297   0.100  1.0000
 10 - 104 -0.170239794 0.20613091 297  -0.826  0.9996
 12 - 16   0.091692745 0.04361807 297   2.102  0.6216
 12 - 20   0.136266428 0.05363752 297   2.541  0.3192
 12 - 26   0.308616325 0.07375087 297   4.185  0.0022
 12 - 32   0.199907724 0.09529330 297   2.098  0.6248
 12 - 44   0.138034505 0.13479517 297   1.024  0.9971
 12 - 52   0.219068800 0.15614450 297   1.403  0.9623
 12 - 64   0.027052521 0.17965340 297   0.151  1.0000
 12 - 78  -0.058577830 0.19570155 297  -0.299  1.0000
 12 - 104 -0.248582297 0.20588322 297  -1.207  0.9881
 16 - 20   0.044573683 0.04361807 297   1.022  0.9971
 16 - 26   0.216923580 0.05993849 297   3.619  0.0177
 16 - 32   0.108214978 0.08092640 297   1.337  0.9735
 16 - 44   0.046341760 0.12249476 297   0.378  1.0000
 16 - 52   0.127376055 0.14603437 297   0.872  0.9993
 16 - 64  -0.064640225 0.17292481 297  -0.374  1.0000
 16 - 78  -0.150270575 0.19217335 297  -0.782  0.9998
 16 - 104 -0.340275043 0.20524946 297  -1.658  0.8855
 20 - 26   0.172349897 0.04807703 297   3.585  0.0199
 20 - 32   0.063641296 0.06670815 297   0.954  0.9984
 20 - 44   0.001768077 0.10925899 297   0.016  1.0000
 20 - 52   0.082802372 0.13479517 297   0.614  1.0000
 20 - 64  -0.109213908 0.16510318 297  -0.661  1.0000
 20 - 78  -0.194844258 0.18784840 297  -1.037  0.9967
 20 - 104 -0.384848725 0.20438557 297  -1.883  0.7690
 26 - 32  -0.108708602 0.04807703 297  -2.261  0.5068
 26 - 44  -0.170581820 0.08813258 297  -1.936  0.7359
 26 - 52  -0.089547525 0.11598318 297  -0.772  0.9998
 26 - 64  -0.281563805 0.15123292 297  -1.862  0.7817
 26 - 78  -0.367194155 0.17965340 297  -2.044  0.6629
 26 - 104 -0.557198623 0.20251567 297  -2.751  0.2069
 32 - 44  -0.061873218 0.06670815 297  -0.928  0.9988
 32 - 52   0.019161077 0.09529330 297   0.201  1.0000
 32 - 64  -0.172855203 0.13479517 297  -1.282  0.9808
 32 - 78  -0.258485554 0.16915399 297  -1.528  0.9317
 32 - 104 -0.448490021 0.19972270 297  -2.246  0.5180
 44 - 52   0.081034295 0.05363752 297   1.511  0.9367
 44 - 64  -0.110981985 0.09529330 297  -1.165  0.9911
 44 - 78  -0.196612335 0.14055285 297  -1.399  0.9631
 44 - 104 -0.386616803 0.19011677 297  -2.034  0.6701
 52 - 64  -0.192016280 0.06670815 297  -2.878  0.1542
 52 - 78  -0.277646630 0.11598318 297  -2.394  0.4135
 52 - 104 -0.467651098 0.17965340 297  -2.603  0.2828
 64 - 78  -0.085630351 0.07375087 297  -1.161  0.9914
 64 - 104 -0.275634818 0.15614450 297  -1.765  0.8354
 78 - 104 -0.190004467 0.11598318 297  -1.638  0.8935

Results are averaged over the levels of: smoking_status 
P value adjustment: tukey method for comparing a family of 12 estimates 

And I see that the backtransformation is gone in the pairwise comparison. I cried another tear

Ideally, I would like to compare the weekly contrasts by the different protein quartile groups, say, compare if contrast week8 - week104 for protein intake group 1 is significantly different thant contrast week8 - week104 for protein intake group 4. This would be equivalent to testing the change score between groups.

Help

Now I'm stuck looking at all this, not being able to extract the estimates, se's, sd's, and confidence intervals that I need(?).

So I was thinking, do I really need to go through all this to answer my question? Or would a simple t-test between change scores for each group suffice?

For some of my outcomes, the protein quartile groups show no effect in the model reduction phase, and therefore gets excluded from the parsimonious statistical model. Should I not reduce the model this far, for the sole purpose of being able to extract estimates (although non-significant) from my statistical model that I can report in my paper?

Since my paper is due in the near future, I need to use the approach that makes most sense to answer my questions efficiently. So should I stick to the emmeans approach (or something similar that works), or should I simply do a two tailed t-test of unadjusted change scores to get my estimates and test sizes?

If I should stick with emmeans, what do I need to do to get the emmeans summary output to reflect the same estimates as the lme summary output. And how do I get the back-transformations to work for the pairs() function? And is there a way to test change scores (post-pre or week104-week8) between the protein groups with the emmeans package as well?

$\endgroup$
1
$\begingroup$

A couple of things regarding emmeans: First, the model you fitted is additive (no interaction terms). This means that there's not much sense in breaking the per-week means down by protein -- especially in doing comparisons. Note that the t statistics and P values for each comparison of weeks are identical with each protein; and that is why. If you do think that the comparisons should be different, then the interaction should be in the model. If you found the interaction was unnecessary in the model, you can just compare weeks marginally.

Second, it helps to understand the order in which things happen. The type = "response" option is an option for the summary() method; it is applied after all the statistics are compiled, just before the summaries are displayed. So with the transformation present, the pairwise comparisons are computed first, and the transformed scale upon which the EMMs lie has no natural counterpart for comparisons among them; that is why the transformation is discarded. If you want differences of back-transformed EMMs, that is possible; the method is to create a new reference grid that is on the transformed scale. This is done using regrid():

 new.emms <- regrid(emmeans.wght_cox)   # transforms everything
 pairs(new.emms)

See the vignette on transformations for details and examples.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.