[R] contrasts with an imbalance in a factor
Connolly, Colm
cgconnolly at ucsd.edu
Fri May 11 21:14:46 CEST 2012
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
More information about the R-help
mailing list