ldply(aovlist, function(x) LSD.test(x, "trt"))
Alan Severini
Dear manipulatr list,
I want to plot bars with errorbars in ggplot2, where the ymin and ymax aesthetics in geom_errorbar will be lower and upper limits of 95% CI calculated as Fisher's LSD (given by the function LSD.test, from package agricolae). So first I need to build a data frame with mean, LCL and UCL. This is similar to Brandon Hurr's thread posted here (http://groups.google.com/group/ggplot2/browse_thread/thread/926079fec134fcce/374b7929457eec2e), since I'm doing an aov function inside dlply(), and then I'm trying to retrieve means, LCL and UCL with ldply(). But instead of performing these functions over an 'aov' or 'lm' object, I need to do it over an 'aovlist' object.
Please, consider the example bellow:
##################################################
calid.isu <- structure(list(rep = c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L,
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L,
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L,
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L,
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), gen = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("B73", "IDS69", "IDS91",
"Mo17", "N209", "R18"), class = "factor"), type = structure(c(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, 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), .Label = c("Dent", "Popcorn"), class = "factor"),
dens = c(3L, 3L, 3L, 9L, 9L, 9L, 3L, 3L, 3L, 9L, 9L, 9L,
3L, 3L, 3L, 9L, 9L, 9L, 3L, 3L, 3L, 9L, 9L, 9L, 3L, 3L, 3L,
9L, 9L, 9L, 3L, 3L, 3L, 9L, 9L, 9L, 3L, 3L, 3L, 9L, 9L, 9L,
3L, 3L, 3L, 9L, 9L, 9L, 3L, 3L, 3L, 9L, 9L, 9L, 3L, 3L, 3L,
9L, 9L, 9L, 3L, 3L, 3L, 9L, 9L, 9L, 3L, 3L, 3L, 9L, 9L, 9L
), trt = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L), .Label = c("Control", "Def17"), class = "factor"),
kw = c(238.0769231, 237.4177778, 224.1810256, 229.8107692,
224.8451282, 228.6366667, 175.2102222, 197.436, 158.1409524,
145.137619, 164.8204762, 144.2605128, 116.3457143, 115.4897778,
112.2661538, 115.6686325, 119.0711111, 112.6721212, 103.1605556,
101.1982222, 97.67285714, 96.57866667, 100.0968889, 94.56888889,
123.6644444, 125.6758974, 125.3346667, 135.3804444, 135.6302222,
132.5013333, 110.8777778, 114.4466667, 115.3431111, 105.8924444,
109.7822222, 104.2968889, 327.79, 328.8592593, 327.37, 262.2291667,
265.9612121, 285.132, 228.1810256, 230.1205128, 272.9977778,
196.8106667, 172.7995556, 173.205641, 216.5764103, 237.9977778,
232.2666667, 182.7433333, 187.4225641, 208.5097778, 141.3138095,
152.408, 136.0119048, 117.8569231, 129.8848889, 120.9033333,
106.4153846, 102.852381, 103.6057143, 97.62153846, 96.84777778,
100.1509524, 92.53055556, 91.49333333, 96.19022222, 82.24761905,
83.12266667, 88.20153846), mc = c(23.25760492, 32.02871633,
34.65023317, 32.45690292, 29.19573376, 35.85711916, 45.62459438,
36.50432267, 45.08864768, 40.89791775, 41.55363895, 50.86208737,
27.73036634, 29.61946087, 28.34098007, 28.11203766, 27.88309525,
30.85823568, 33.9165003, 30.66898858, 27.74011773, 34.70989204,
30.95391683, 25.81645751, 30.23371769, 29.13251504, 29.18093537,
27.83517939, 24.00401265, 28.1040084, 31.66388188, 31.37946578,
27.65534832, 32.83527251, 29.33132122, 32.97439847, 32.52239547,
35.92882075, 34.68844361, 39.73504983, 34.2269868, 34.15030242,
49.54189863, 45.5438863, 42.04615753, 40.00304928, 51.32046204,
50.618281, 48.76762857, 45.08466849, 45.03301055, 45.84655552,
47.73670599, 38.3534471, 58.7864382, 56.65173149, 56.45083267,
61.27120314, 59.79169157, 58.4880375, 22.9867415, 32.31126237,
26.619675, 28.64651744, 29.36804256, 29.44483249, 16.84879612,
24.97032454, 29.69108211, 33.58592578, 30.3770339, 32.51187632
), pkw = c(242.2768126, 228.680255, 239.4746315, 262.039447,
263.2004701, 226.8675081, 243.6634693, 242.0392536, 248.6912731,
228.4701575, 237.3306151, 246.0736283, 152.3000756, 158.0675412,
155.6719293, 151.9645413, 148.2571533, 153.6609524, 158.4086807,
158.0675412, 155.6719293, 149.5316815, 148.2571533, 151.9725864,
162.5885414, 160.9788268, 158.8064109, 152.3719954, 172.9538174,
167.7573828, 161.4546129, 169.800479, 161.2746472, 154.6772647,
172.9538174, 167.7573828, 389.4870901, 384.3005886, 335.4383425,
309.4714107, 270.21271, 293.7798344, 396.8964478, 337.8219796,
325.1478363, 282.7629727, 290.6703963, 256.2895962, 317.3198801,
301.936277, 303.5155856, 269.1139456, 278.8664954, 303.8879498,
308.6469766, 303.461522, 306.004106, 299.5171102, 280.7100043,
264.5428868, 152.2062599, 154.5055613, 148.5175863, 153.3532503,
155.2284616, 143.7431684, 152.2062599, 154.5055613, 148.5175863,
153.3532503, 142.0665971, 149.7421524), pgrf = c(282.960257,
254.1734275, 240.6155626, 197.680563, 139.5904623, 139.6596155,
47.62081779, 73.00704218, -96.22175785, 24.76535785, 45.40737471,
22.0900658, 138.3327921, 106.9787992, 135.2926594, 112.2698487,
89.24703803, 54.16606432, 2.094338429, -15.30589435, 11.76090324,
-32.94815199, -56.60310487, -72.11302741, 110.3469882, 121.447115,
128.8139126, 77.51340409, 48.38849394, 74.60649737, -31.59691349,
-36.9670944, -51.98043774, -41.71178562, -61.63155052, -52.78620601,
210.5014205, 159.9687629, 196.3622316, 89.84514132, 107.0525139,
124.4276131, 17.57997959, 20.29328792, -21.31351317, 21.97909617,
1.199363404, -42.23004236, 285.9688313, 246.1021303, 275.44742,
91.13503225, 131.7465335, 164.0287927, 19.44279222, 81.73469842,
50.81670903, -22.71307468, 12.70481589, -14.45791459, 104.6696135,
167.8091461, 166.4906192, 54.34342927, 44.84734735, 111.1693072,
11.23854043, 101.6383664, -0.042240981, 1.986861097, -19.9889307,
-10.89982554), ssr = c(0.377970946, 0.378668462, 0.337479512,
0.406651309, 0.400836675, 0.357683824, 0.072959677, 0.136022815,
-0.171056761, 0.054567916, 0.096407366, 0.046990272, 0.413301108,
0.270175614, 0.368271643, 0.331202034, 0.294132424, 0.181534742,
-0.000835446, -0.038283974, 0.025715447, -0.106157921, -0.196447722,
-0.381516188, 0.27132605, 0.337224988, 0.329679261, 0.222336839,
0.156485734, 0.23099572, -0.081055715, -0.106325213, -0.147177335,
-0.164378341, -0.41218537, -0.23971544, 0.463786116, 0.413196239,
0.487843653, 0.265630026, 0.417197671, 0.357140499, 0.040526958,
0.051952848, 0.063802296, 0.075893275, 0.009891788, -0.072509733,
0.361686174, 0.329997471, 0.416584006, 0.228905458, 0.295672446,
0.356176786, 0.037527093, 0.099007957, 0.074649034, -0.413948461,
0.025148999, -0.081669963, 0.2084008, 0.186590857, 0.198032412,
0.034684933, 0.0915229, 0.157155298, 0.00613107, 0.12973367,
0.006071002, -0.016209485, -0.088115066, -0.026633937), pp = c(8.226666667,
8.78, 8.893333333, 7.4, 7.464285714, 7.546666667, 9.913333333,
10.32666667, 10.6, 8.607142857, 7.946666667, 7.728571429,
10.72142857, 10.66666667, 10.27333333, 10.12, 9.966666667,
9.62, 10.16153846, 10.15333333, 9.88, 9.32, 8.873333333,
9.218181818, 10.57333333, 10.79230769, 10.27333333, 10.15333333,
9.586666667, 10, 10.22, 10.42, 9.846666667, 8.921428571,
8.469230769, 9.291666667, 9.028571429, 9.213333333, 9.378571429,
8.3, 8.257142857, 7.88, 9.593333333, 10.04, 10.27857143,
9.246666667, 8.88, 8.4, 9.62, 9.906666667, 9.628571429, 9.107142857,
9.371428571, 8.635714286, 11.84666667, 11.21333333, 11.68,
10.47692308, 10.12666667, 10.72307692, 13.09567308, 12.00982143,
12.67008929, 13.31730769, 12.48461538, 12.93375, 11.34427083,
11.85708333, 12.24333333, 11.25535714, 11.48208333, 11.00714286
), pw = c(19.63068923, 20.82280489, 20.07296359, 16.95479385,
16.87866769, 17.52835444, 17.37838533, 20.36367867, 16.52403476,
12.5863419, 13.25804714, 11.16846769, 12.58136, 12.32354133,
11.46933282, 11.66742286, 11.86551289, 10.82495818, 10.49355333,
10.30237156, 9.67978619, 8.990371111, 8.876934667, 8.609257576,
13.01994167, 13.54327897, 12.86824178, 13.73020267, 12.99338578,
13.23400533, 11.34608311, 11.9667119, 11.35500844, 9.585918095,
9.226869744, 10.25411889, 29.92234111, 31.5053437, 31.112615,
21.497105, 22.27312933, 22.31708667, 21.96421846, 23.19897897,
28.36303833, 18.21953156, 15.332492, 14.85881556, 20.77974308,
23.62982278, 22.37470333, 16.41339636, 17.59638256, 18.13831238,
16.76006238, 17.04711467, 15.81661905, 12.49370556, 13.107072,
13.27119444, 13.96985375, 12.42366813, 13.11851827, 13.03034157,
12.22666639, 12.97252756, 10.49532167, 10.83901522, 11.80635636,
9.255307708, 9.796418417, 9.783658109), yield = c(157.2776927,
139.9024787, 149.9553147, 97.027102, 64.74249733, 73.41981867,
92.72405467, 92.10446733, 83.13818133, 49.31235733, 60.32463933,
47.39602467, 32.34655533, 38.983946, 36.464944, 34.23528367,
32.00562333, 28.61189733, 31.24101867, 31.96166, 28.66819,
26.37701467, 27.145134, 16.02664, 46.27985667, 44.87895267,
44.69626133, 39.381704, 35.13821533, 35.79286467, 39.31776867,
34.92770733, 35.21303333, 25.04913, 18.26496267, 19.083308,
121.5008153, 109.358096, 120.5726333, 73.584356, 67.600308,
87.27384933, 69.62566267, 75.09868533, 78.319708, 47.30491267,
39.76569667, 37.59053133, 148.16205, 151.7672793, 162.0824393,
57.00391467, 69.64438533, 82.30209067, 69.28333067, 94.60482933,
72.357718, 32.14987067, 43.233688, 34.55512267, 47.07733333,
94.67579333, 69.49366667, 37.88815333, 40.05880667, 52.17333333,
51.3877, 58.25456, 50.93452, 33.1968, 29.88622, 23.02785333
)), .Names = c("rep", "gen", "type", "dens", "trt", "kw",
"mc", "pkw", "pgrf", "ssr", "pp", "pw", "yield"), row.names = c(NA,
72L), class = "data.frame")
library(reshape)
calid.isu <- melt(calid.isu, id.vars = c("rep",
"gen",
"type",
"dens",
"trt"))
for (i in 1:5) { calid.isu[, i] <- as.factor(calid.isu[, i])
}
library(plyr)
calid.isu.aov <- dlply(calid.isu,
"variable",
function(x) aov(value ~ rep +
dens +
trt +
type +
type/gen +
type:dens +
type/gen:dens +
type:trt +
type/gen:trt +
dens:trt +
type:dens:trt +
type/gen:dens:trt +
Error(rep/dens/trt/gen),
data = x))
library(agricolae)
lsd <- ldply(calid.isu.aov, function(x) LSD.test(x,
"dens",
alpha = 0.05,
p.adj = "none"))
##################################################
The last command gives error, since LSD.test can't work with an 'aovlist' object. Looking inside the structure of calid.isu.aov, it looks that one can access each of the stratums of the aovlist (each of one is an aov) by a code like the following:
calid.isu.aov$kw[[3]]
So I presume that if one need a data frame of means, LCL and UCL for the 3rd stratum (where the error term is 'rep:dens' and the treatment is 'dens') the code could be:
lsd <- ldply(calid.isu.aov, function(x) LSD.test(x[[3]],
"dens",
alpha = 0.05,
p.adj = "none"))
But it doesn't work either. I'm sorry because this question may be naive, as a result of a misunderstanding of the basic structure of a list object, or how plyr works, but can someone give me some clue about how to solve this problem?
Thanks in advance,
Alan
Brandon Hurr
Brandon Hurr
Luciano Selzer
--
You received this message because you are subscribed to the Google Groups "manipulatr" group.
To post to this group, send email to manip...@googlegroups.com.
To unsubscribe from this group, send email to manipulatr+...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/manipulatr?hl=en.
Luciano Selzer
A <- y$model
There is no model in this type of anova. So it returns NULL and breaks everything up. Which means that a function in a package I mantain is also broken in this situation...
Alan Severini
Thank you very much for your help. So I will contact the agricolae
package administrator. I wonder is there is another function to do
multiple comparisons that could work with an aovlist.
Cheers,
Alan
El 03/08/2011 05:00 p.m., Luciano Selzer escribi�:
> Follow up, I've found the bug.
>
>
>
>
> A<- y$model
>
>
>
> There is no model in this type of anova. So it returns NULL and breaks
> everything up. Which means that a function in a package I mantain is
> also broken in this situation...
> Luciano
>
>
> 2011/8/3 Luciano Selzer <luciano...@gmail.com
> <mailto:luciano...@gmail.com>>
>
> Indeed, it is related to the error term because there are like
> diferent aovs inside that object, each one whit it's own error term.
> I thougth that doing LSD.test(kw.aov[[3]], trt="dens") would work.
> But it returns an empty response. So I think that maybe LSD.test
> doesn't work or we are doing something wrong. You could try to fix
> it yourself or contact the package maintainer.
>
> HTH
>
> Luciano
>
>
>
> 2011/8/3 Brandon Hurr <brando...@gmail.com
> <mailto:brando...@gmail.com>>
> <alanse...@hotmail.com <mailto:alanse...@hotmail.com>>
> wrote:
>
> aov(value ~ rep +
> dens +
> trt +
> type +
> type/gen +
> type:dens +
>
> type/gen:dens +
> type:trt +
>
> type/gen:trt +
> dens:trt +
>
> type:dens:trt +
>
> type/gen:dens:trt +
>
> Error(rep/dens/trt/gen),
> data = x))
>
>
>
> --
> You received this message because you are subscribed to the
> Google Groups "manipulatr" group.
> To post to this group, send email to manip...@googlegroups.com
> <mailto:manip...@googlegroups.com>.
> To unsubscribe from this group, send email to
> manipulatr+...@googlegroups.com
> <mailto:manipulatr%2Bunsu...@googlegroups.com>.
Dennis Murphy
I had a similar experience to Brandon and Luciano, but I have concerns
with the model you're using. According to your error term, you are
modeling a five-level hierarchy: rep/dens/trt/gen + residual.
According to this specification, you have treatment nested within dens
and gen nested within rep * dens * trt combinations. The table
generated by the line below tells a different story:
with(calid.isu, ftable(gen ~ rep + type + dens + trt))
The eight replicates correspond to the number of levels of 'variable',
used to separate by-groups.
This table shows that rep is a true replication factor: the within rep
design is the same in all three reps. So far, so good.
Within rep, we find that levels B73, Mo17 and N209 only occur in
combination with type Dent, and levels IDS69, IDS91 and R18 occur in
combination with type Popcorn. Therefore, gen is nested within Type.
At each level of gen nested within type, dens and trt form a 2 x 2
factorial structure. Therefore, this is a split-plot design with an
RCB structure and nesting within the levels of type at the whole-plot
level and a 2 x 2 CRD at the gen level. The corresponding ANOVA should
reflect this structure. Taking a data subset corresponding to the
variable 'kw' and naming it tst, I tried a model in the lme4 package
that seemed more or less reasonable (not having plotted the data, I
can't be sure, but the structure looks plausible). OTOH, I have no
sense of context of the data, using only the table above as a guide to
its structure.
> mm <- lmer(value ~ rep + type * dens * trt + (0 + type|gen), data = tst)
> summary(mm)
Linear mixed model fit by REML
Formula: value ~ rep + type * dens * trt + (0 + type | gen)
Data: tst
AIC BIC logLik deviance REMLdev
551.6 583.4 -261.8 575.8 523.6
Random effects:
Groups Name Variance Std.Dev. Corr
gen typeDent 1842.47 42.924
typePopcorn 141.46 11.894 0.314
Residual 147.70 12.153
Number of obs: 72, groups: gen, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 261.673 25.193 10.387
rep2 3.232 3.508 0.921
rep3 1.929 3.508 0.550
typePopcorn -148.765 26.346 -5.647
dens9 -32.805 5.729 -5.726
trtDef17 -75.413 5.729 -13.163
typePopcorn:dens9 34.349 8.102 4.240
typePopcorn:trtDef17 63.331 8.102 7.817
dens9:trtDef17 -3.433 8.102 -0.424
typePopcorn:dens9:trtDef17 -4.569 11.458 -0.399
<Correlation matrix of fixed effects estimates snipped>
You can use the multcomp package on objects of class lmer (like mm),
which provide multiple comparison procedures that are far more
accurate than LSD. The statistical problem with looking at only the
level associated with Type is that the variance estimate is likely to
be wrong. One has to be very careful to get the correct variance
estimates of fixed effects and contrasts in hierarchical models
because variance propagates up the hierarchy - the expected mean
squares show how the individual variance components are to be
combined.
Bottom line suggestion:
(1) Use the lme4 package to fit the models;
(2) Use the multcomp package to produce the multiple comparisons;
(3) Write one function to fit the model and the multiple comparisons;
extract output from the multiple comparisons run and put it into a
data frame. It's easier if this function is written outside the ddply()
call. It should take a (generic) data frame as input and a data
frame as output.
(4) Call ddply() with variable as the grouping factor and the function
in (3) as the third argument.
HTH,
Dennis
Alan Severini
Dear Dennis,
Thanks for your deep analysis of this problem and specially for the bottom line guidelines. I'm afraid that given my level of ignorance (e.g. what does "(0 + type | gen)" mean in the lmer model?) I cannot answer most of your concerns. But you are right in that:
- this is a split-plot design
- gen is nested within type
I see that the summary of an lmer model is very different from an aov one. So, in order to extract the some kind of 95 % CI by using the package multcomp on lmer models, I must try hard to learn about the lme4 package.
Best regards,
Alan
----------------------------------------
> Date: Wed, 3 Aug 2011 14:36:51 -0700
> Subject: Re: ldply(aovlist, function(x) LSD.test(x, "trt"))
> From: djm...@gmail.com
> To: alanse...@hotmail.com
> CC: manip...@googlegroups.com