atorus-research / CDISC_pilot_replication

A modern replication of the CDISC pilot table outputs, using the PHUSE Test Data Factory ADaM data, using the R programming language
https://www.atorusresearch.com
MIT License
46 stars 22 forks source link

MMRM analysis #6

Open slophaven opened 3 years ago

slophaven commented 3 years ago

I have been looking a bit into the MMRM analysis in the report for "CDISCPILOT01 – Initial Case Study of the CDISC SDTM/ADaM Pilot Project". More specifically I looked into the generation of table 14-3.11 and the SAS output supporting table 14-3.11. The title of the table is "... change from baseline to week 24". It turns out that the numbers presented are not for the change to week 24, it is the average change to week 8/16/24. The numbers presented are created using this piece of SAS code (after having applied the following filter: EFFFL == "Y" & PARAMCD == 'ACTOT' & ANL01FL == 'Y' & DTYPE != 'LOCF' & AVISITN > 0):

proc mixed noclprint data=temp01 ic ; class trtpcd_f aweekc sitegr1 usubjid; model chg=trtpcd_f aweekc sitegr1 trtpcd_faweekc base baseaweekc/s ddfm=kr; repeated aweekc/subject=usubjid type=un; lsmeans trtpcd_f/ diff cl; ods output diffs=temp02; ods output LSMeans=temp03; run;

What should be done instead is this (note the difference in the lsmeans statement):

proc mixed noclprint data=temp01 ic ; class trtpcd_f aweekc sitegr1 usubjid; model chg=trtpcd_f aweekc sitegr1 trtpcd_faweekc base baseaweekc/s ddfm=kr; repeated aweekc/subject=usubjid type=un; lsmeans trtpcd_f*aweekc/ diff cl; ods output diffs=temp02; ods output LSMeans=temp03; run;

In R I did the following:

First approach

test1<-lmer(CHG ~ TRTPCD_F + SITEGR1 + AWEEKC + TRTPCD_F:AWEEKC + BASE + BASE:AWEEKC + (AVISITN | USUBJID),data=adas)

Second approach using gls

test1<-gls(CHG ~ TRTPCD_F + SITEGR1 + AWEEKC + TRTPCD_F:AWEEKC + BASE + BASE:AWEEKC, corr=corSymm(form= ~1 | USUBJID),weights=varIdent(form= ~1 | AVISITN),data=adas)

Third approach: lmer with modification as suggested by Daniel Sabanes Bove

test1<-lmer(CHG ~ TRTPCD_F + SITEGR1 + AWEEKC + TRTPCD_F:AWEEKC + BASE + BASE:AWEEKC + (0 + AWEEKC | USUBJID),data=adas, control=lmerControl(check.nobs.vs.nRE="ignore"))

LSmeans per group

test2<-emmeans::lsmeans(test1, ~TRTPCD_F:AWEEKC, lmer.df='kenward-roger')

LSmeans for differences

test3<-emmeans::contrast(test2, method="pairwise", adjust=NULL)

95% CI

test4 <- confint(test3)

The results are as follows:

In SAS:

PBO: 2.3291 [0.9680; 3.6903] High dose: 1.5009 [-0.1475; 3.1494] Low dose: 1.7352 [0.2247; 3.2457] Differences: High vs. PBO: 0.8282 [-1.2856; 2.9420] p=0.4403 Low vs. PBO: 0.5939 [-1.4136; 2.6014] p=0.5600

First approach: with lmer:

Pbo Week 24 2.32 [0.962; 3.68] Xan_Hi Week 24 1.49 [-0.154; 3.14] Xan_Lo Week 24 1.77 [0.263; 3.28] Differences: Xan_Hi Week 24 - Pbo Week 24 0.8270 [-1.283; 2.9367] p=0.4404 Xan_Lo Week 24 - Pbo Week 24 0.5459 [-1.459; 2.5510] p=0.5919

Second approach: with gls (using approximative Satterthwaite):

Pbo Week 24 2.331 [0.972; 3.69] Xan_Hi Week 24 1.502 [-0.140; 3.14] Xan_Lo Week 24 1.735 [0.229; 3.24] Differences: Xan_Hi Week 24 - Pbo Week 24 0.8294 [-1.278; 2.9371] p=0.4384 Xan_Lo Week 24 - Pbo Week 24 0.5959 [-1.407; 2.5991] P=0.5578

Third approach: lmer with modification as suggested by Daniel Sabanes Bove

Pbo Week 24 2.329 [0.968; 3.69] Xan_Hi Week 24 1.501 [-0.147; 3.15] Xan_Lo Week 24 1.735 [0.225; 3.25] Differences: Xan_Hi Week 24 - Pbo Week 24 0.8282 [-1.285; 2.9415] p=0.4402 Xan_Lo Week 24 - Pbo Week 24 0.5939 [-1.413; 2.6010] p=0.5599