Dr. Winsemius,
Really quick, a BACI is a Before-After-Control-Impact approach. I have a long
time series of sea cucumber density estimates, which are taken at the same
location(s) through time. Some are in areas Impacted by sea otters and some
are in areas Not Impacted by sea otters (two levels). Each estimate is also
coded with a "B" if it is Before the time sea otters showed up at the impacted
sites or "A" if after (Impact and Control sites are both coded with BA). BACI
analyses suggest an impact if the ANOVA interaction term (BA*Otter) is
significant; i.e., changes in sea cucumber density from before to after depend
on whether sea otters are present. I log transformed the response to help
normalize the data, as it has many zeros. While shapiro-wilks does not suggest
normality, it is a very large data set and it is "approximately" normal based
on graphical examination. Again, the data is unbalanced, as there are many more
estimates at the control sites and before period.
With that said... I would like to perform the following pairwise comparisons;
B:With Otters v. A:With Otters, B:Without Otters v. A:Without Otters
I am performing other ANOVAs with different data and no interaction, where I
would like to perform multiple pairwise comparisons between the fivelevels of a
single factor. I used the code I provided previously and still received error
messages. If I can get this BACI interaction problem figured out, I should
manage to adjust it to other models. I recently came accross Dunnett's
Modified Tukey Kramer (DTK.test) and it appears to address the same issue of
unbalanced data and has a very straightforward script (although I am not sure
it lends itself to interactions?). Is this test an appropriate substitute for
the glht method?
You wrote, "When your factors are both binary, the effect estimates fit nicely
into a 2 x 2 table and the consideration of the single contrast added by the
interaction is fairly simple" wool=='A'
wool=='B' tension=='L' 3.7179 -0.4356+3.7179
tension=='M' -0.6012+3.7179 0.6281+3.7179
This output seems to be what I am looking for; assuming that if the value range
for a comparison includes zero, then there is no significant difference? Where
did those values come from?I hope this helps clear up my problem. If you have
concerns about pitfalls with this approach, then I would love to hear them and
I can research them outside of this thread. This is part of a masters thesis
and needs to be sound.
Thank you very much for your time.
Sean
Date: Thu, 14 Jun 2012 10:26:58 -0700
From: ml-node+s789695n4633417...@n4.nabble.com
To: seanlars...@hotmail.com
Subject: Re: Tukey Kramer with ANOVA (glm)
On Jun 13, 2012, at 7:36 PM, Alaska_Man wrote:
> Hello,
>
> I am performing a BACI analysis with ANOVA using the following glm:
I admit I had no idea what a "BACI analysis" might be. Looking it up
it appears to be a cross-over design and my statistical betters have
sternly warned me about this regression briar patch in the past. I'm
especially suspicious of the lack of any statements about the balance
in the sampling in your presentation. (And for that matter the
extremely sketchy statement of design.)
>
> fit1<-glm(log(Cucs_m+1)~(BA*Otter)+BA+Otter+ID+Primary, data=b1)
I'm guessing you do not understand that BA*Otter in an R formula
expands to BA + Otter + BA:Otter
> The summary(aov(fit1)) shows significance in the interaction;
> however, now I
> would like to determine what combinations of BA and Otter are
> significantly
> different (each factor has two levels). ID and PRIMARY substrates are
> categorical and included in the model to help explain some of the
> variation
> in the data. The data is unbalanced so I plan on using Tukey Kramer
> post
> hoc analysis. Here is how my data is laid out, it is a fairly
> substantial
> data set:
Editing done on original (although it proved unrevealing.)
> Subdistrict T Year Cucs_m Primary Persistence Otter
> Fishing BA ID
> 109-41,42 9 2010 0.00 sil 3 1
> 1 A 109-41,42
> 109-41,42 13 2010 2.75 rck 3 1
> 1 A 109-41,42
> 109-41,42 16 2010 2.00 rck 3 0
> 1 A 109-41,42
> 109-41,42 18 2010 8.25 rck 3 0
> 0 B 109-41,42
>
> I am assuming this is an appropriate pairwise comparison analysis
> and I
> cannot get the code to work with my data.
What does it mean to be doing "pairwise comparisons" on two-level
factor variables?)
> I am *unclear how to code it to
> work with the interaction*; however, even when I attempt to use it
> only for
> a single factor, it does not work (see below).
>
> x<-aov(glm(Cucs_m~as.factor(BA),data=cuc))
> glht(x, linfct=mcp(BA="Tukey"))
> ....................................
> Error in mcp2matrix(model, linfct = linfct) :
> Variable(s) âBAâ have been specified in âlinfctâ but cannot be
> found in
> âmodelâ!
I suspect the glht() function is looking for 'as.factor(BA)` in the
model matrix and not finding it. If BA is not already a factor, then
it would make sense to do:
cuc$BA <- factor(cuc$BA)
.... before any analysis. Notice that you get a warning that
performing contrasts in the presence of interactions is something to
be warned about. If you do not know what you are doing here (and your
proposed analysis hints at that possibility), I may have set a trap
for you by solving a syntactic problem but not solving a conceptual
problem.
> mod <- glm(log(breaks) ~ wool*tension, data=subset(warpbreaks,
tension %in% c("L","M")))
> glht(mod, linfct=mcp(tension="Tukey"))
General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Linear Hypotheses:
Estimate
M - L == 0 -0.6012
Warning message:
In mcp2matrix(model, linfct = linfct) :
covariate interactions found -- default contrast might be
inappropriate
---------------
Looking at the mod-object you see that the "Estimate" above is
actually NOT what you had interest in. You were presumably more
interested in the contrast woolB:tensionM whose coefficient was 0.6281.
----
Coefficients:
(Intercept) woolB tensionM woolB:tensionM
3.7179 -0.4356 -0.6012 0.6281
----------
I would have instead done something like this:
> mod <- glm(log(breaks) ~ wool*tension, data=subset(warpbreaks,
tension %in% c("L","M")))
> mod2 <- glm(log(breaks) ~ wool+tension, data=subset(warpbreaks,
tension %in% c("L","M")))
> anova(mod,mod2)
Analysis of Deviance Table
Model 1: log(breaks) ~ wool * tension
Model 2: log(breaks) ~ wool + tension
Resid. Df Resid. Dev Df Deviance
1 32 4.6235
2 33 5.5113 -1 -0.88777
Now I can say that the addition of an interaction term resulted in a
non-significant improvement in model fit at least when measured on the
log(breaks) scale. (Note: This is quite a different result than one
sees on the untransformed scale where the interaction is highly
significant.) When your factors are both binary, the effect estimates
fit nicely into a 2 x 2 table and the consideration of the single
contrast added by the interaction is fairly simple.
wool=='A' wool=='B'
tension=='L' 3.7179 -0.4356+3.7179
tension=='M' -0.6012+3.7179 0.6281+3.7179
> Can anyone off suggestions on potential problems with my approach
> and/or
> script issues?
Why was the log transformation being done? Is the desired outcome a
statement about ratios?
--
David Winsemius, MD
West Hartford, CT
______________________________________________
[hidden email] 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.
If you reply to this email, your message will be added to the
discussion below:
http://r.789695.n4.nabble.com/Tukey-Kramer-with-ANOVA-glm-tp4633314p4633417.html
To unsubscribe from Tukey Kramer with ANOVA (glm), click here.
NAML
--
View this message in context:
http://r.789695.n4.nabble.com/Tukey-Kramer-with-ANOVA-glm-tp4633314p4633435.html
Sent from the R help mailing list archive at Nabble.com.
[[alternative HTML version deleted]]
______________________________________________
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.
