Hi everybody,
I have an experiment examining risky choice behavior where two groups of
subjects were unevenly divided across two different MRI scanners while they
performed a task. Each subject's data was recorded once and only once on a
particular scanner. The table describing the distribution of subjects across
the scanner (3TE and 3TW) and groups is below.
3TE 3TW
Group1 10 9
Group2 3 18
I'm trying to perform linear mixed effects analysis to determine the effect of
group (Group1 and Group2) and task.
I can run the lme without problems, however, when I then go to perform a number
of contrasts using the contrast package it barfs because of the unequal number
of rows in the scanner factor.
My question is if there a way for me to specify contrasts (with a package like
contrast) in a way that can deal with the fact that one of the factors in my
model is unevenly balanced?
The code below illustrates the problem I'm having.
Any help would be most gratefully appreciated.
-------------------------------------------->8------------------------------------------------------------
rm(list=ls())
library(nlme)
library(contrast)
my.model=structure(list(
subject = structure(c(1L, 2L, 3L, 10L, 11L, 12L,
15L, 18L, 20L, 27L, 31L, 34L, 36L, 35L, 33L, 38L, 39L, 37L, 40L,
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, 38L, 39L, 37L, 40L,
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, 38L, 39L, 37L, 40L,
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, 38L, 39L, 37L, 40L,
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, 38L, 39L, 37L, 40L,
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", "T0244",
"T0245", "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, 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, 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, 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, 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, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L),
.Label = c("Group1", "Group2"), class = "factor"),
task = 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,
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, 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),
.Label = c("20outcome",
"40outcome", "80outcome", "neg40outcome", "neg80outcome"), class =
"factor"),
scanner = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
1L, 1L, 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, 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,
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, 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, 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.0444161854684353,
-0.114592373371124, 0.0769545212388039, -0.268393993377686,
0.0290516708046198, -0.0482161603868008, 0.123185224831104,
0.0947600975632668, -0.157943934202194, 0.129291623830795,
-0.0649581402540207, 0.0572351925075054, -0.163955569267273,
0.00085166294593364, 0.0297060012817383, -0.0351637043058872,
-0.0439938083291054, -0.0113321952521801, 0.0744024813175201,
0.0953737944364548, -0.00487061077728868, -0.0259095914661884,
0.0361779294908047, -0.27682438492775, 0.00985431391745806,
0.130222201347351, -0.0171694327145815, -0.108230642974377,
-0.024471502751112, -0.0748439505696297, 0.00296992133371532,
-0.113702893257141, -8.38618052512174e-06, -0.134054526686668,
0.0275573581457138, -0.0549034662544727, 0.0316105224192142,
0.0187845807522535, 0.0888368785381317, 0.0392464064061642,
0.00777457281947136, -0.0786111205816269, 0.034312792122364,
-0.431242674589157, 0.0550877340137959, -0.0284439660608768,
0.14971898496151, 0.18965645134449, -0.249666750431061, -0.114646390080452,
0.120740003883839, 0.0412555411458015, 0, 0.0195636544376612,
0.0360089540481567, 0.141455665230751, -0.106815293431282,
-0.060357641428709, 0, 0.032408569008112, 0.0176739227026701,
0, 0.0896040201187134, -0.340650588274002, 0.298701554536819,
0.0824474170804024, 0.0104535883292556, -0.100377805531025,
0, 0.0286602880805731, 0.055004108697176, -0.149470701813698,
-0.0354152917861938, 0.0355890430510044, 0.0701439455151558,
0, -0.0408579632639885, -0.0490981005132198, 0, 0.118443757295609,
0, 0.0385815761983395, -0.0136463027447462, 0.995342969894409,
0.0263149011880159, 0.0687721744179726, -0.0559106767177582,
0, -0.248157367110252, -0.763005673885345, -0.198830619454384,
0.108319185674191, 0.148703306913376, -0.0386029928922653,
0.0344023033976555, 0.225042834877968, 0.0763957425951958,
-0.0267908964306116, 0, 0.0849404707551003, -0.22040219604969,
0, 0.198203936219215, 0.156402811408043, -0.170966222882271,
0, 0.113942489027977, -0.0061397822573781, 0.177627995610237,
0.00290966127067804, 0, -0.121305607259274, 0.0242178849875927,
0.288999557495117, 0.151954352855682, 0, -0.00922275241464376,
0.113115973770618, 0, 0.183086037635803, 0.0219641979783773,
0.086994431912899, 0.00716547947376966, -0.156700059771538,
-0.0983942002058029, 0.0657249838113785, 0.16248020529747,
0.065130889415741, -0.177818566560745, 0.0923013314604759,
-0.338354974985123, 0.168124347925186, 0.126458197832108,
-0.140643388032913, -0.0198399741202593, 0.082248218357563,
0.1715097874403, -0.149834543466568, 0, 0.177407532930374,
0.0246455036103725, 0, -0.0490625686943531, -0.0734507888555527,
0.289667189121246, 0.0715472921729088, 0.412098109722137,
-0.0883435308933258, 0.0285825077444315, -0.211081445217133,
0.101080782711506, -0.0608764216303825, -0.0586727559566498,
-0.188042506575584, -0.0544104352593422, -0.0368366353213787,
0.133244857192039, -0.150556549429893, 0, 0.205774813890457,
0, -0.0538171827793121, -0.0553216300904751, 1.29600214958191,
0.0555772967636585, -0.181610956788063, 0.262428969144821,
0, -0.107578098773956, 0.196385011076927, 0.0106251891702414,
0.218737199902534, 0, -0.198070287704468, -0.142257764935493,
0.277058124542236, 0.033035684376955, -0.104592069983482,
0, -0.0898942425847054, 0.0699265524744987, 0, -0.0476390048861504,
0.181289047002792, 0.146841675043106, 0, 0.0353749208152294,
-0.0641172826290131, 0.328668653964996, -0.0261893272399902,
-0.0553865134716034, -0.00882275588810444, 0.0129395183175802,
0, 0.128713548183441, 0, -0.203418999910355, 0.0572037324309349,
0, -0.0326705500483513)),
.Names = c("subject", "group",
"task", "scanner", "fmri"), class = "data.frame", row.names = c(NA,
-200L))
## this is what I really want to be able to do, but without contrast
## complaining about the imbalance in the number of rows
cat("### With scanner\n")
fixedFormula=as.formula("fmri ~ group * task + scanner")
randomFormula = as.formula("random = ~ 1 | subject")
mylme = lme(fixed=fixedFormula, random=randomFormula, data=my.model)
## now look at the risky choices (40outcome and 80outcome) versus the
## safe choices (20outcome)
con=contrast(
mylme,
a=list(group="Group1", task=c("40outcome", "80outcome"),
scanner=levels(my.model$scanner)),
b=list(group="Group1", task="20outcome",
scanner=levels(my.model$scanner)))
print(con)
## the code below wont execute becasue the contrast call above will
## error out
### not what I want but just to show that it works when scanner is not
### includede in the model
cat("### Without scanner\n")
fixedFormula=as.formula("fmri ~ group * task")
randomFormula = as.formula("random = ~ 1 | subject")
mylme = lme(fixed=fixedFormula, random=randomFormula, data=my.model)
con=contrast(
mylme,
a=list(group="Group1", task=c("40outcome", "80outcome")),
b=list(group="Group1", task="20outcome"))
print(con)
Regards,
--
Colm G. Connolly, Ph. D.
Dept of Psychiatry
University of California, San Diego
9500 Gilman Dr. #0738
La Jolla, CA 92093-0738
Tel: +1-858-246-0684
______________________________________________
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.
