Hi,
I have a problem with running a script in parallel with the foreach function
for computing linear models at every voxel of a set of 3D images. When I run
it using %do%, it runs fine and I get the correct output but when I run it
with %dopar% then my main output is not correct. The problem inside the
foreach loop seems to be in the function called "lmer" which is part of the
lmerTest package. It is used to compute linear models and the issue is that
lmer returns the p-values which are different for each test when using %do%
but they are all the same when using %dopar%, which is obviously incorrect.
Here is the foreach loop:
#To try with the first 5 voxels only, so that it works.
resultTest<-foreach (index=1:5, .combine=rbind, .packages='lmerTest',
.verbose=TRUE) %dopar% {
i1 <- inmask[index,1]
i2 <- inmask[index,2]
i3 <- inmask[index,3]
y <- ydat[i1,i2,i3,]
# create the datafile with current voxel, with the voxelvalues and
variables of
# interest
long <- data.frame(cbind(xdat,y))
# do the model you want to do on the data
fm1 <- lmer(y ~ 1 + time + (1|id), long)
# extract what you care about and save in your created array, with the
index you are
# currently running. Type fm1 to see your output, str(fm1) to see your
slots, so you
# can find where the data you want are (which place in the order, which
can
depend on
# your model).
p=f...@t.pval[2]
p.time[i1,i2,i3] <- p
}
And here is the output p-values when running with %do%:
result.1 0.02976868
result.2 0.01706419
result.3 0.01509690
result.4 0.02078050
result.5 0.03950854
Output using %dopar% (I ran this one right after the %do%, you can see that
all the p-values here are the same as the last p-value ran with %do%)
result.1 0.03950854
result.2 0.03950854
result.3 0.03950854
result.4 0.03950854
result.5 0.03950854
Also, here is the lmer function (from the lmerTest package):
lmer <-function(formula, data, family = NULL, REML = TRUE,
control = list(), start = NULL, verbose = FALSE, doFit = TRUE,
subset, weights, na.action, offset, contrasts = NULL,
model = TRUE, x = TRUE, ...)
{
mc<-match.call()
mc[[1]] <- quote(lme4::lmer)
model<-eval.parent(mc)
model<-as(model,"merLmerTest")
#tryCatch( { result = glm( y~x , family = binomial( link = "logit" )
) } , error = function(e) { print("test") } )
t.pval <- tryCatch( {totalAnovaRandLsmeans(model=model,
ddf="Satterthwaite", isTtest=TRUE)$ttest$tpvalue}, error = function(e) {
NULL })
if(!is.null(t.pval))
{
mo...@t.pval <-t.pval
}
else
{
model<-as(model,"mer")
}
return(model)
}
What I do not understand is that all the other data output by lmer is fine
and different for each voxel except for this p-value.
So if anyone has a hint on how I could solve this issue, I would greatly
appreciate.
Thank you in advance.
Best regards,
Alexis.
--
View this message in context:
http://r.789695.n4.nabble.com/Issue-with-lmerTest-and-dopar-tp4676831.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
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.
