s-news
[Top] [All Lists]

Re: lme and aov for for fractional design with a random effect

To: Olivier Renaud <Olivier.Renaud@pse.unige.ch>
Subject: Re: lme and aov for for fractional design with a random effect
From: Spencer Graves <spencer.graves@pdf.com>
Date: Thu, 19 Feb 2004 07:37:54 -0800
Cc: Snews <s-news@wubios.wustl.edu>
In-reply-to: <6.0.2.0.0.20040219160237.026c3820@mail.pse.unige.ch>
References: <6.0.2.0.0.20040219160237.026c3820@mail.pse.unige.ch>
User-agent: Mozilla/5.0 (Windows; U; Windows NT 5.0; en-US; rv:1.4) Gecko/20030624 Netscape/7.1 (ax)
It is my impression that short, clear questions are more likely to be answered than longer, more complicated questions. It helps, at least for me, if the question contains a toy example in the form of S code that I can copy out of an email and test for myself. That helps me quickly test a suggestion before replying with a stupid answer than doesn't work. For "aov" have you consulted Venables and Ripley (2002) Modern Applied Statistics with S (Springer)? For "lme", have you consulted Pinheiro had Bates (2000) Mixed-Effects Models in S and S-Plus (Springer)? Also, have you consulted the R archives at "www.r-project.org" -> search -> "R site search"? I'm sorry I can't answer your question, but I hope this will help. spencer graves

Olivier Renaud wrote:

I asked the following question 10 days ago, and received no answer. If you 
don't have time for a lengthy answer, a pointer to a reference would be 
appreciated.
Olivier

Hi,
I'm not exactly sure how to call the design that we use, but it may contain the 
keywords fractional design, repeated measure, BIB, MLM, and maybe more! Here is 
the setting: We have three types of route (a1, a2, a3) and three conditions 
(b1, b2, b3). Since each subject (S) can only do each route once, each subject 
was asked to do only 3 experiments:
for 1/6 of subjects, it was the 3 combinations (a1,b1), (a2,b2), (a3,b3)
for 1/6 of subjects, it was the 3 combinations (a1,b1), (a2,b3), (a3,b2)
for 1/6 of subjects, it was the 3 combinations (a1,b2), (a2,b1), (a3,b3)
for 1/6 of subjects, it was the 3 combinations (a1,b2), (a2,b3), (a3,b1)
for 1/6 of subjects, it was the 3 combinations (a1,b3), (a2,b1), (a3,b2)
for 1/6 of subjects, it was the 3 combinations (a1,b3), (a2,b2), (a3,b1)

I guess the formulas for aov and lme are (on simulated data)
summary (aov (resp ~ A+B+Error(S/(A+B)), data=esfr5))
anova (lme (resp ~ A+B, random=~1|S, data=esfr5), type="marginal")
# See Output and Data below

with a * instead of the + if one includes the interactions. I have the 
following questions:

1) Why is the error in aov divided only in two pieces. If the design were complete (each 
subject having done 9 experiments), it would be divided in four? Moreover, why is the 2nd 
Error term called "A %in% S"?

2) Is the interaction A:B estimable or is it confounded with something? As you 
noticed, it appears twice in the tests. What do these two tests mean? This is a 
simulated example, but what if they are very different?

3) Up to now, with _real_ data, lme always find singularity ('Problem in 
.C("mixed_EM",: Singularity in backsolve, while calling subroutine mixed_EM' ). 
Are they some parameters in lmeControl that can be tuned to solve this ? (aov has no 
problem with these data, see output of one example of real data and the data below)

4) Are the adjusted models different between lme and aov (especially when 
including interaction) ?

5) Can we analyse these data with SPSS ?

Olivier







#Output of simulated data

summary(aov(resp~A+B+Error(S/(A+B)),data=esfr5))
Error: S Df Sum of Sq Mean Sq F Value Pr(F) Residuals 29 3.114277 0.1073889 Error: A %in% S Df Sum of Sq Mean Sq F Value Pr(F) A 2 0.071461 0.0357305 0.458006 0.6348918
       B  2  0.206622 0.1033109 1.324277 0.2741954
Residuals 56 4.368734 0.0780131
summary(aov(resp~A*B+Error(S/(A*B)),data=esfr5))
Error: S Df Sum of Sq Mean Sq F Value Pr(F) A:B 4 0.491633 0.1229084 1.171608 0.347148 Residuals 25 2.622643 0.1049057 Error: A %in% S Df Sum of Sq Mean Sq F Value Pr(F) A 2 0.071461 0.0357305 0.463168 0.6318575
       B  2  0.206622 0.1033109 1.339203 0.2709411
     A:B  4  0.357267 0.0893167 1.157798 0.3401417
Residuals 52 4.011467 0.0771436
anova(lme(resp~A+B, random=~1|S, data=esfr5), type="marginal")
numDF denDF F-value p-value (Intercept) 1 56 190.1546 <.0001
         A     2    56   0.4580  0.6349
         B     2    56   1.3243  0.2742
anova(lme(resp~A*B, random=~1|S, data=esfr5), type="marginal")
numDF denDF F-value p-value (Intercept) 1 52 190.4814 <.0001
         A     2    52   0.4597  0.6340
         B     2    52   1.3292  0.2735
       A:B     4    52   1.0644  0.3835


# Output of real data:
summary(aov(Reponse~Lien+Repetition+Error(Sujet/(Lien+Repetition)),data=SPenSLR1))
Error: Sujet Df Sum of Sq Mean Sq F Value Pr(F) Residuals 23 4984396 216712.9 Error: Lien %in% Sujet Df Sum of Sq Mean Sq F Value Pr(F) Lien 2 62062.6 31031.30 2.917880 0.0645736
Repetition  2   35207.5 17603.76 1.655285 0.2027117
Residuals 44 467934.7 10634.88
anova(lme( fixed=Reponse~Lien+Repetition, data=SPenSLR1, random= ~1| Sujet), 
type="marginal")
Problem in .C("mixed_EM",: Singularity in backsolve, while calling subroutine mixed_EM Use traceback() to see the call stack



# Simulated and real data:

"esfr5" <- structure(.Data = list(c("s1", "s1", "s1", "s2", "s2", "s2", "s3", "s3", "s3", "s4", "s4", "s4", "s5", "s5", "s5", "s6", "s6", "s6", "s7", "s7", "s7",
       "s8", "s8", "s8", "s9", "s9", "s9", "s10", "s10", "s10", "s11", "s11", "s11", "s12", "s12", "s12", 
"s13", "s13", "s13",
       "s14", "s14", "s14", "s15", "s15", "s15", "s16", "s16", "s16", "s17", "s17", "s17", "s18", "s18", 
"s18", "s19", "s19",
       "s19", "s20", "s20", "s20", "s21", "s21", "s21", "s22", "s22", "s22", "s23", "s23", "s23", "s24", 
"s24", "s24", "s25",
       "s25", "s25", "s26", "s26", "s26", "s27", "s27", "s27", "s28", "s28", "s28", "s29", "s29", "s29", 
"s30", "s30", "s30")
, "value" = c("a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", 
"a3", "a1", "a2", "a3",
       "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", 
"a3", "a1", "a2",
       "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", 
"a2", "a3", "a1",
       "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3", 
"a1", "a2", "a3",
       "a1", "a2", "a3", "a1", "a2", "a3", "a1", "a2", "a3")
, "value" = structure(.Data = c(1, 2, 3, 1, 3, 2, 2, 1, 3, 2, 3, 1, 3, 1, 2, 3, 
2, 1, 1, 2, 3, 1, 3, 2, 2, 1, 3, 2, 3, 1, 3, 1, 2, 3, 2, 1, 1, 2, 3, 1, 3,
       2, 2, 1, 3, 2, 3, 1, 3, 1, 2, 3, 2, 1, 1, 2, 3, 1, 3, 2, 2, 1, 3, 2, 3, 
1, 3, 1, 2, 3, 2, 1, 1, 2, 3, 1, 3, 2, 2, 1,
       3, 2, 3, 1, 3, 1, 2, 3, 2, 1)
, levels = c("b1", "b2", "b3")
, class = "factor"
)
, c(0.024646461009979248, 0.82596908044070005, 0.39935583062469959, 
0.42794170277193189, 0.79330063052475452, 0.41740808542817831,
0.42689276579767466, 0.97277549048885703, 0.86306705698370934, 0.20230660401284695, 0.33371580624952912, 0.61999470554292202, 0.74001671373844147, 0.07387995719909668, 0.30246197292581201, 0.35847964184358716, 0.46099079307168722, 0.16043389076367021, 0.58081981865689158, 0.22667973302304745, 0.16983740823343396, 0.33375463215634227, 0.41230997582897544, 0.35489231906831264, 0.1146096084266901, 0.66199403256177902, 0.086105858441442251, 0.92259258404374123, 0.033672890160232782, 0.97718893829733133, 0.54872243478894234, 0.063163440208882093, 0.18481955723837018, 0.69088406767696142, 0.00048982584848999977, 0.9480914156883955, 0.18231043219566345, 0.006615795660763979, 0.67624678136780858, 0.31703903619199991, 0.38895514793694019, 0.21916917618364096, 0.80299417721107602, 0.61142651597037911, 0.60505218058824539, 0.75811595609411597, 0.59731168439611793, 0.90336279338225722, 0.66943846317008138, 0.86255600396543741, 0.58649370679631829, 0.64013873925432563, 0.35643625492230058, 0.1302819293923676, 0.76756644900888205, 0.39251976786181331, 0.19385471474379301, 0.01424878928810358, 0.21007986413314939, 0.059990861918777227, 0.7896508164703846, 0.25749274902045727, 0.77087316988036036, 0.50042055826634169, 0.96990229655057192, 0.95752570079639554, 0.68353205313906074, 0.7501282668672502, 0.41614685719832778, 0.24187999684363604, 0.68989087454974651, 0.80914638098329306, 0.44907701155170798, 0.83931859489530325, 0.93146510748192668, 0.44827951770275831, 0.66668012505397201, 0.82564923400059342, 0.090691552963107824, 0.61973605211824179, 0.74060451937839389, 0.033790021203458309, 0.62395105883479118, 0.56522526312619448, 0.32841544877737761, 0.67579050175845623, 0.029812175780534744, 0.0094875660724937916, 0.4210558719933033, 0.067989643197506666)
)
, row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", 
"19", "20", "21", "22",
       "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", 
"40", "41", "42",
       "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", 
"60", "61", "62",
       "63", "64", "65", "66", "67", "68", "69", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79", 
"80", "81", "82",
       "83", "84", "85", "86", "87", "88", "89", "90")
, names = c("S", "A", "B", "resp")
, class = "data.frame"
)



"SPenSLR1" <- structure(.Data = list(c(1164.5, 968., 931.25, 1263.25, 1091., 1189.5, 802., 655.5, 689.75, 1067.5, 1178.75, 1370.75, 909.5, 878.5, 1070., 649.25,
       594.25, 580.75, 769., 776.75, 688., 1489.75, 1361.25, 1428.25, 718.5, 
608.25, 612.5, 892.25, 1509.75, 1372.75, 1477.5,
       1284., 1425.5, 1014.25, 981.5, 1218.5, 798.75, 840.25, 1022.75, 966.75, 
1049., 1218., 1212.5, 1027.5, 1149.5, 798.75,
       938.25, 919., 1034.5, 947.5, 1028.5, 997., 900.25, 1120.25, 611.5, 688., 
621.75, 680.25, 722., 685.5, 822.75, 914.,
       973.75, 912.75, 969., 992.75, 1602.25, 1596.5, 1716.25, 1200.25, 1281., 
1275.)
, structure(.Data = c(7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 
12, 12, 13, 13, 13, 14, 14, 14, 1, 1, 1, 15, 15, 15, 16, 16, 16, 17,
       17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 
23, 23, 2, 2, 2, 24, 24, 24, 25, 25, 25, 26,
       26, 26, 27, 27, 27, 28, 28, 28)
, levels = c("2", "3", "4", "6", "7", "9", "10", "11", "12", "13", "15", "17", "18", "19", "20", "21", "22", "23", 
"24", "25", "26", "27",
       "29", "31", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", 
"49", "50", "51",
       "52", "53", "55", "59", "60", "61", "62", "64", "65")
, class = "factor"
)
, structure(.Data = c(1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 
2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2,
       3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 
3, 1, 2, 3, 1, 2, 3)
, levels = c("1", "2", "3", "4")
, class = "factor"
)
, structure(.Data = c(1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 3, 2, 1, 3, 2, 1, 
3, 2, 1, 3, 2, 2, 1, 3, 2, 1, 3, 2, 1, 3, 2, 1, 3, 2, 3, 1, 2, 3,
       1, 2, 3, 1, 2, 3, 1, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 2, 1, 3, 2, 
1, 3, 2, 1, 3, 2, 1)
, levels = c("1", "2", "3", "4")
, class = "factor"
)
)
, names = c("Reponse", "Sujet", "Lien", "Repetition")
, class = "data.frame"
, row.names = c("1", "6", "11", "17", "22", "27", "33", "38", "43", "49", "54", "59", "65", "71", "74", "81", "87", 
"90", "97", "103", "106",
       "113", "119", "122", "130", "133", "139", "146", "149", "155", "162", "165", "171", "178", "181", 
"187", "194", "199",
       "201", "210", "215", "217", "226", "231", "233", "242", "247", "249", "259", "261", "266", "275", 
"277", "282", "291",
       "293", "298", "307", "309", "314", "323", "326", "329", "339", "342", "345", "355", "358", "361", 
"371", "374", "377")
)
--------------------------------------------------------------------
This message was distributed by s-news@lists.biostat.wustl.edu.  To
unsubscribe send e-mail to s-news-request@lists.biostat.wustl.edu with
the BODY of the message:  unsubscribe s-news


<Prev in Thread] Current Thread [Next in Thread>