Hi everybody,
I've written a function to run an LME model on data derived from functional
magnetic resonance images. When I run the function with contrasts included I get
the following error
Error in eval(expr, envir, enclos) : object 'inModelFormula' not found
I think it has something do do with the way contrast evaluates arguments, but
I've got no idea how to fix it. The code attached below demonstrates the
problem. Could anyone advise on how to fix/eliminate the problem?
------------------------------------------------------------------------->8----------------------------------------------------------------------
rm(list=ls())
library(nlme)
library(contrast)
model = structure(list(subject = structure(c(1L, 2L, 3L, 10L, 11L, 12L,
15L, 18L, 20L, 27L, 31L, 34L, 36L, 35L, 33L, 37L, 40L,
39L, 41L,
38L, 42L, 5L, 6L, 7L, 8L, 13L, 22L, 14L, 21L, 19L, 25L,
23L,
24L, 16L, 28L, 26L, 30L, 9L, 29L, 32L, 17L, 4L, 1L, 2L,
3L, 10L,
11L, 12L, 15L, 18L, 20L, 27L, 31L, 34L, 36L, 35L, 33L,
37L, 40L,
39L, 41L, 38L, 42L, 5L, 6L, 7L, 8L, 13L, 22L, 14L, 21L,
19L,
25L, 23L, 24L, 16L, 28L, 26L, 30L, 9L, 29L, 32L, 17L,
4L, 1L,
2L, 3L, 10L, 11L, 12L, 15L, 18L, 20L, 27L, 31L, 34L,
36L, 35L,
33L, 37L, 40L, 39L, 41L, 38L, 42L, 5L, 6L, 7L, 8L, 13L,
22L,
14L, 21L, 19L, 25L, 23L, 24L, 16L, 28L, 26L, 30L, 9L,
29L, 32L,
17L, 4L, 1L, 2L, 3L, 10L, 11L, 12L, 15L, 18L, 20L, 27L,
31L,
34L, 36L, 35L, 33L, 37L, 40L, 39L, 41L, 38L, 42L, 5L,
6L, 7L,
8L, 13L, 22L, 14L, 21L, 19L, 25L, 23L, 24L, 16L, 28L,
26L, 30L,
9L, 29L, 32L, 17L, 4L, 1L, 2L, 3L, 10L, 11L, 12L, 15L,
18L, 20L,
27L, 31L, 34L, 36L, 35L, 33L, 37L, 40L, 39L, 41L, 38L,
42L, 5L,
6L, 7L, 8L, 13L, 22L, 14L, 21L, 19L, 25L, 23L, 24L,
16L, 28L,
26L, 30L, 9L, 29L, 32L, 17L, 4L), .Label =
c("T0004", "T0005",
"T0009",
"T0020", "T0026", "T0030", "T0049",
"T0050", "T0072",
"T0078",
"T0079", "T0080", "T0083", "T0085",
"T0094", "T0104",
"T0105",
"T0111", "T0112", "T0113", "T0114",
"T0115", "T0119",
"T0131",
"T0136", "T0141", "T0143", "T0150",
"T0152", "T0162",
"T0167",
"T0196", "T0198", "T0216", "T0217",
"T0218", "T0219",
"T0244",
"T0245", "T0246", "T0253", "T0262"),
class = "factor"),
group = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L), .Label = c("A", "B"), class =
"factor"),
brik = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L), .Label = c("20r", "40r",
"80r", "neg40r",
"neg80r"), class = "factor"),
scanner = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L), .Label = c("3TE", "3TW"), class
= "factor"),
fmri = c(0.0274474211037159, -0.293868184089661, -0.0320024862885475,
-0.293503433465958, 0.143126413226128, 0.030569875612855,
0.292777240276337, 0.0842535793781281, 0.370734214782715,
-0.205730959773064, -0.159094661474228, -0.0162228047847748,
0.0576154813170433, 0.0974738970398903, 0.0753356739878654,
0.383649528026581, 0.152102336287498, -0.0200147740542889,
0.563074350357056, 0.0299493446946144, 0.224797010421753,
-0.137106865644455, 0.113709755241871, -0.0426492169499397,
0.125936061143875, -0.0884832739830017, -0.184091582894325,
-0.0396079421043396, 0.137704193592072, -0.16732420027256,
0.323838174343109, 0.0527379289269447, 0.16394031047821,
0.0882217213511467, -0.086088553071022, 0.263536393642426,
-0.0624078810214996, 0.501975774765015, 0.293505430221558,
0.0178689565509558, -0.0358020290732384, -0.0210242234170437,
-0.0215112138539553, -0.190718099474907, -0.184049472212791,
0.129670202732086, 0.127691984176636, 0.103092163801193,
0.028329448774457, -0.0824636742472649, 0.239311337471008,
0.052397083491087, -0.243118211627007, -0.0604248456656933,
0, 0.0283053312450647, 0.0447858348488808, 0.153856858611107,
0.161434963345528, 0.133259132504463, 0.345422327518463,
0.179489523172379, 0, 0.0857120081782341, -0.227173984050751,
0, 0.108787253499031, -0.0804603472352028, 0.303540110588074,
0.0349666140973568, -0.0674800649285316, -0.218899860978127,
0, -0.0647929981350899, -0.0807700902223587, 0.189496427774429,
-0.14031058549881, 0.102083601057529, 0.0112999016419053,
0, 0.146465063095093, -0.0535526983439922, 0, -0.0516969263553619,
0, -0.0553484782576561, -0.0141696771606803, -0.234366849064827,
-0.10783002525568, -0.476115167140961, -0.159414410591125,
0, -0.034053809940815, -0.325537651777267, -0.0498526841402054,
-0.455243200063705, -0.138358011841774, -0.00459326291456819,
0.0563298426568508, 0, 0.0871630012989044, 0.0294053312391043,
-0.145188868045807, 0.185977503657341, 0, 0.0104816136881709,
-0.422605484724045, 0, 0.0189609378576279, 0.0802518352866173,
-0.332018971443176, 0, 0.1592066437006, 0.015329628251493,
0.464863181114197, -0.0143416309729218, 0, 0.0231552720069885,
-0.0517721027135849, 0.246574491262436, 0.040972862392664,
0, 0.0851080566644669, 0.100515216588974, 0, 0.0195400360971689,
-0.0863960310816765, -0.12854839861393, -0.247136428952217,
0.120662644505501, 0.276547312736511, 0.203144460916519,
0.223532348871231, 0.0516968779265881, 0.538363635540009,
0.212853536009789, 0.230947688221931, -0.479386448860168,
-0.155570849776268, -0.0206140652298927, -0.286562293767929,
-0.136779427528381, -0.183586418628693, -0.0615491755306721,
-0.193362891674042, 0.350995421409607, 0, 0.193982109427452,
0.337202131748199, 0, -0.337008684873581, -0.0792389065027237,
-0.894769728183746, 0.0251937098801136, 0.0576306283473969,
-0.0896522775292397, -0.769122302532196, -0.0115023702383041,
0.0801776126027107, 0.0963655114173889, 0.0338999405503273,
0.147212103009224, -0.101382903754711, 0.548872232437134,
0.146211251616478, -0.0617087669670582, 0, 0.181398138403893,
0, -0.196751445531845, -0.179890334606171, 0.0863914489746094,
0.192345842719078, 0.0343081317842007, -0.219587966799736,
0, 0.0491224192082882, 0.315233290195465, 0.188764750957489,
-0.340905874967575, 0, 0.154585316777229, -0.13356277346611,
0, -0.0309638362377882, 0.241891011595726, 0.108909144997597,
0.0359761491417885, 0, 0.168808549642563, 0.0375387519598007,
0, 0.0292902421206236, 0.0512374006211758, 0.185121148824692,
0, 0.330981016159058, -0.221254393458366, 0.624710321426392,
0.328806430101395, -0.138702377676964, -0.045757669955492,
-0.188088029623032, 0, 0.0608050674200058, 0, 0.42381739616394,
0.115315444767475, 0, 0.201624438166618)), .Names = c("subject",
"group",
"brik", "scanner", "fmri"), row.names = c(NA,
-210L), class = "data.frame")
runLme = function (inData, inZ, inNumberOfOutputBriks, inContrasts,
inAnovaIndices, inModel, inModelFormula, inRandomFormula) {
cat("There will be " , inNumberOfOutputBriks, " stats
briks\n")
outStats <- vector(mode="numeric", length=inNumberOfOutputBriks)
##cat (inData, "\n")
if ( ! all(inData == 0 ) ) {
##try(
##if( inherits(
print(inModelFormula)
print(inRandomFormula)
mylme <- lme(fixed=inModelFormula, random=inRandomFormula, data =
inModel)
print(mylme)
##,
## silent=FALSE),
##"try-error") ) {
##temp <- 0
## cat (paste("Error on slice", inZ, "\n"))
##} else {
temp <- as.vector(unlist(anova(mylme, type="marginal")))
##}
if(length(temp) > 1) {
numberOfContrasts=length(inContrasts)
outStats[1:(inNumberOfOutputBriks-(numberOfContrasts*2))] =
temp[c(inAnovaIndices)]
if (numberOfContrasts > 0 ) {
contrastsStartAt=(inNumberOfOutputBriks-(numberOfContrasts*2))+1
cat("Contrasts start at ", contrastsStartAt, "\n")
for (i in seq(1, numberOfContrasts, by=1)) {
print( inContrasts[[i]])
con1=contrast(mylme, a =inContrasts[[i]]$a, b =inContrasts[[i]]$b,
type="average")
if (length(con1) > 1) {
outStats[contrastsStartAt:contrastsStartAt+1] <- c(con1$Contrast,
con1$testStat)
}
contrastsStartAt=contrastsStartAt+2
}
}
}
} else {
cat("Got all zeros :-(\n")
}
return(outStats)
}
scannerLevels=c("3TE", "3TW")
prePlannedContrasts=list(
"response" = list(
list(
name="aRiskyVaSafe",
a=list("group"="A",
"brik"=c("40r"), "scanner"=scannerLevels),
b=list("group"="A",
"brik"=c("20r"),
"scanner"=scannerLevels)),
list(
name="bRiskyVbSafe",
a=list("group"="B",
"brik"=c("40r"), "scanner"=scannerLevels),
b=list("group"="B",
"brik"=c("20r"),
"scanner"=scannerLevels)),
list(
name="aRiskyVbRisky",
a=list("group"="A",
"brik"=c("40r", "80r"),
"scanner"=scannerLevels),
b=list("group"="B",
"brik"=c("40r", "80r"),
"scanner"=scannerLevels)),
list(
name="aSafeVbSafe",
a=list("group"="A",
"brik"=c("20r"), "scanner"=scannerLevels),
b=list("group"="B",
"brik"=c("20r"), "scanner"=scannerLevels))
)
)
## the formulae used for the fixed effects (modelFormula) and random effects
(randomFormula)
modelFormula = as.formula("fmri ~ group * brik + scanner")
randomFormula = as.formula("random = ~ 1 | subject")
anovaIndices=c(1L, 6L, 11L, 16L, 2L, 7L, 12L, 17L, 3L, 8L, 13L, 18L, 4L, 9L,
14L, 19L, 5L, 10L, 15L, 20L)
numberOfOutputBriks=28
## formal parameters
## (inData, inZ, inNumberOfOutputBriks, inContrasts,
inAnovaIndices, inModel, inModelFormula, inRandomFormula)
a= runLme(model$fmri, 21, numberOfOutputBriks,
prePlannedContrasts[["response"]], anovaIndices, model,
modelFormula, randomFormula)
> version
_
platform x86_64-redhat-linux-gnu
arch x86_64
os linux-gnu
system x86_64, linux-gnu
status
major 2
minor 14.1
year 2011
month 12
day 22
svn rev 57956
language R
version.string R version 2.14.1 (2011-12-22)
> packageVersion("contrast")
[1] ?0.16?> packageVersion("nlme")
[1] ?3.1.102?
Regards,
--
Colm G. Connolly, Ph. D.
Dept of Psychiatry
University of California, San Diego
9500 Gilman Dr. #0738
La Jolla, CA 92093-0738