I fit the following model in R using the lme4 package, which is a linear mixed effects regression of outcome "WT_abs_change" (absolute weight change) on various factors, including "ns(WT_dtime, df=3)*regimen" (interaction between time, using a natural/restricted cubic spline, and drug regimen, which is categorical with 5 categories):
fm1_adj<-lmer(WT_abs_change ~ ns(WT_dtime, df=3)*regimen + WT_baseline + ns(age_med, df=3) + female + black + MSM + ns(reg1yr, df=3) + ns(CD4_baseline, df=3) + (1|CFAR_PID), data=WT_naive_dataset) summary(fm1_adj)
The model output follows:
summary(fm1_adj) Linear mixed model fit by REML ['lmerMod'] Formula: WT_abs_change ~ ns(WT_dtime, df = 3) * regimen + WT_baseline +
ns(age_med, df = 3) + female + black + MSM + ns(reg1yr, df = 3) +
ns(CD4_baseline, df = 3) + (1 | CFAR_PID) Data: WT_naive_dataset
REML criterion at convergence: 71605.9
Scaled residuals: Min 1Q Median 3Q Max -16.6026 -0.4862 0.0077 0.4769 5.8467
Random effects:
Groups Name Variance Std.Dev.
CFAR_PID (Intercept) 79.84 8.935
Residual 81.74 9.041
Number of obs: 9559, groups: CFAR_PID, 1174
Fixed effects: Estimate Std. Error t value (Intercept) 10.659383 2.781723 3.832 ns(WT_dtime, df = 3)1 15.618293 1.440984 10.839 ns(WT_dtime, df = 3)2 17.506589 1.685360 10.387 ns(WT_dtime, df = 3)3 13.294438 1.126701 11.799 regimenEVG 0.358854 1.562922 0.230 regimenNNRTI -0.876973 1.830318 -0.479 regimenPI -1.064086 1.755134 -0.606 regimenRAL -0.985957 2.187114 -0.451 WT_baseline -0.009162 0.006922 -1.324 ns(age_med, df = 3)1 1.991276 1.361413 1.463 ns(age_med, df = 3)2 0.297468 3.127315 0.095 ns(age_med, df = 3)3 1.097583 2.426688 0.452 female -0.087786 0.983952 -0.089 black 0.341828 0.605586 0.564 MSM 0.559518 0.789013 0.709 ns(reg1yr, df = 3)1 2.067972 1.320659 1.566 ns(reg1yr, df = 3)2 -3.173422 2.460025 -1.290 ns(reg1yr, df = 3)3 1.123731 1.874192 0.600 ns(CD4_baseline, df = 3)1 -8.578638 1.964845 -4.366 ns(CD4_baseline, df = 3)2 -20.852911 3.082316 -6.765 ns(CD4_baseline, df = 3)3 -7.753365 4.906481 -1.580 ns(WT_dtime, df = 3)1:regimenEVG -9.270812 1.965526 -4.717 ns(WT_dtime, df = 3)2:regimenEVG -14.157056 2.306043 -6.139 ns(WT_dtime, df = 3)3:regimenEVG -11.326906 1.528430 -7.411 ns(WT_dtime, df = 3)1:regimenNNRTI -7.658012 1.697332 -4.512 ns(WT_dtime, df = 3)2:regimenNNRTI -7.420794 1.973527 -3.760 ns(WT_dtime, df = 3)3:regimenNNRTI -6.678979 1.317237 -5.070 ns(WT_dtime, df = 3)1:regimenPI -5.167451 1.622453 -3.185 ns(WT_dtime, df = 3)2:regimenPI -0.628741 1.899374 -0.331 ns(WT_dtime, df = 3)3:regimenPI -5.158171 1.258856 -4.098 ns(WT_dtime, df = 3)1:regimenRAL -5.068909 2.344669 -2.162 ns(WT_dtime, df = 3)2:regimenRAL -0.727831 2.863420 -0.254 ns(WT_dtime, df = 3)3:regimenRAL -7.413167 1.832161 -4.046
I can then also use "visreg" to get lovely plots of outcome over time, by category of "regimen" (using the interaction including splines):
visreg(fm1_adj, "WT_dtime", by="regimen", overlay=TRUE, band=TRUE, partial=FALSE, ylab="Weight (lbs.)", xlab="Months after ART-Start", ylim=c(150,190))
My question is: is there a shortcut or simple way to do a postestimation test (either LRT or Wald test) for differences between categories of my category-by-spline interaction, incorporating the spline terms appropriately (ideally using the model output, but potentially using the visreg predicted results)? Any advice would be appreciated; thank you!