Hi,
I have been running successfully Westfall procedures to adjust p-values
in multiple comparisons.
However, the last model I tried consumes too much memory (>40 Gb).
Is there a way to reduce memory requirements? Or a workaround to run
multiple comparison tests with Westfall procedures?
Mind that I don't fully understand what's going on on the Westfall
procedure implemented. So sorry in advance if this is a stupid question.
Here is the script I am running:
require("nlme")
require("multcomp")
dt_wide <- data.frame(
Subj=1:44,
Status=factor(rep(c("D","S","B","D","S","B"), c(7,7,6,8,8,8))),
Treatment=factor(rep(c("St","Un"), c(20,24))),
day1=rnorm(44,4,0.5),
day2=rnorm(44,4,0.5),
day3=rnorm(44,4,0.5))
dt_wide[,4:6] <- apply(dt_wide[,4:6], 2, function(x) {x[sample(1:44,
sample(0:3))] <- NA;x})
dt_long <- reshape(dt_wide, idvar="Subj", times=names(dt_wide)[4:6],
timevar="Time",
varying=list(names(dt_wide)[4:6]), v.names="Values", direction="long")
dt_long$Int <- with(dt_long, interaction(Status, Time, Treatment))
fm <- lme(Values ~ Int, data=dt_long, random= ~1 |Subj, na.action="na.omit")
K<-rbind("1; St: D-S "=c(0,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0),
"1; St: B-D "=c(0,-1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0),
"1; St: B-S "=c(0, 0,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0),
"1; Un: S-D "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 0,
0, 0, 0),
"1; Un: B-D "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 0, 1, 0, 0,
0, 0, 0),
"1; Un: B-S "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0,
0, 0, 0),
"2; St: S-D "=c(0, 0, 0, 0,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0),
"2; St: B-D "=c(0, 0, 0, 0,-1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0),
"2; St: B-S "=c(0, 0, 0, 0, 0,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0),
"2; Un: S-D "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 1,
0, 0, 0),
"2; Un: B-D "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 0,
1, 0, 0),
"2; Un: B-S "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,
1, 0, 0),
"3; St: S-D "=c(0, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0),
"3; St: B-D "=c(0, 0, 0, 0, 0, 0, 0,-1, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0),
"3; St: B-S "=c(0, 0, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 0, 0, 0,
0, 0, 0),
"3; Un: S-D "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0,-1, 1),
"3; Un: B-D
"=c(0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-2,-1),
"3; Un: B-S
"=c(0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-2),
"D; 1: Un-St"=c(0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0),
"S; 1: Un-St"=c(0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0),
"B; 1: Un-St"=c(0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0),
"D; 2: Un-St"=c(0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0),
"S; 2: Un-St"=c(0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0),
"B; 2: Un-St"=c(0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0),
"D; 3: Un-St"=c(0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0),
"S; 3: Un-St"=c(0, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0,
0, 0, 1),
"B; 3:
Un-St"=c(0,-1,-1,-1,-1,-1,-1,-1,-1,-2,-1,-1,-1,-1,-1,-1,-1,-1),
"D; St: 2-1 "=c(0,-1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0),
"D; Un: 2-1 "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 1, 0,
0, 0, 0),
"S; St: 2-1 "=c(0, 0,-1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0),
"S; Un: 2-1 "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0, 1,
0, 0, 0),
"B; St: 2-1 "=c(0, 0, 0,-1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0),
"B; Un: 2-1 "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 0, 0,
1, 0, 0),
"D; St: 3-2 "=c(0, 0, 0, 0,-1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0),
"D; Un: 3-2 "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 0,
0, 1, 0),
"S; St: 3-2 "=c(0, 0, 0, 0, 0,-1, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0),
"S; Un: 3-2 "=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,
0, 0, 1),
"B; St: 3-2 "=c(0, 0, 0, 0, 0, 0,-1, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0),
"B; Un: 3-2
"=c(0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-2,-1,-1))
colnames(K) <- names(fixef(fm))
mc <- glht(fm, linfct=K)
summary(mc,test=adjusted(type="Westfall"))
Thanks,
Joao
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
