I am analyzing an experimental data set. The data consists of a paired vector of treatment type and a binomial outcome:
Can I use an ANOVA (one-way or factorial?) for this or do I have to use another test? Can anyone recommend a suitable test?
$\begingroup$
$\endgroup$
1
-
1$\begingroup$ Can you clarify if by "paired" you mean each subject experiences each treatment and has an outcome for each treatment? If the answer to that is "No," then Gregg H's answer works. If the answer is "Yes," however, there's a different test. $\endgroup$– AlexisCommented Apr 12, 2018 at 15:35
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
9
This is many months late, but an ANOVA followed by t-tests if you care for multiple comparisons would be appropriate. This may go against conventional wisdom. However, when we are modeling a response variable using regression techniques, we usually care about how well we can model the mean of the variable and how reliable our statistical inference is. The extent to which our ability to infer about relations when the response is at its mean should help us decide the adequacy of a particular model.
As to the first question, regardless of what regression model you choose, logistic, probit, ANOVA, the predicted means of the response on the probability scale will be the exact same values since your single predictor is a grouping variable. So all models will yield identical fit to the response variable in terms of prediction. This changes if you add covariates in an ANCOVA type situation. Here's some sample R code to show predicted probabilities are the same:
set.seed(12345)
# create three category grouping variable
x <- sample(1:3, 100, TRUE)
y <- rep(NA, 100)
# grouping 1 average probability = .2, group 2 = .4, group 3 = .6
y[x == 1] <- rbinom(length(y[x == 1]), 1, .2)
y[x == 2] <- rbinom(length(y[x == 2]), 1, .4)
y[x == 3] <- rbinom(length(y[x == 3]), 1, .6)
x <- factor(x) # make x categorical
# Obtain predicted probabilities
unique(fitted(glm(y ~ x))) # anova
[1] 0.5263158 0.6000000 0.2187500
unique(fitted(glm(y ~ x, binomial))) # logistic
[1] 0.5263158 0.6000000 0.2187500
unique(fitted(glm(y ~ x, binomial(link = "probit")))) # probit
[1] 0.5263158 0.6000000 0.2187500
The next question is about statistical inference. With a binary outcome, your errors and residuals if you check them will neither be normally distributed nor will they have constant variance, so you violate some of the classical assumptions. In practice though, it does not matter. There is a 1972 paper by Glass, Peckham and Sanders that talks about this after a review of the literature. It is also easy to check this by simulation. Your p-values behave like they do when you meet classical assumptions. Moreover, power is not affected relative to using a logistic regression model. It is useful to know that there is a weighted least squares alternative to ANOVA that is technically correct, but in this situation, it does not matter.
Here's a simple simulation to check whether ANOVA maintains the nominal error rate when the null is true. It's a three group design, balanced, with total sample size of 50.
set.seed(12345)
p <- replicate(500, {
x <- factor(sample(1:3, 50, TRUE))
y <- rbinom(50, 1, .5)
# x and y are uncorrelated
# conduct ANOVA and save p values
summary(aov(y ~ x))[[1]]$`Pr(>F)`[1]
})
# what is the rate of false positives?
# it should be close to 5% if the null is true with alpha at .05
mean(p < .05)
[1] 0.042
There are many situations where problems will arise should a linear model be applied to a binary response variable, but ANOVA and t-tests are immune to these problems.
Paper by Glass, Peckham and Sanders: http://journals.sagepub.com/doi/10.3102/00346543042003237
answered Jul 9, 2018 at 12:02
-
$\begingroup$ A few comments. 1) The 5% rate you mentioned is not power. It is the type-I error rate (false positive). The way you conducted the simulation, you aren't able to calculate the power because your simulation has no cases of true positives. 2) Your simulation hasn't addressed what is may be the most salient question, do the two approaches (anova and logistic regression) come to the same conclusion? From the way you did the simulation, the p-values from the anovas track very closely to the p-values from logistic regression, although I'm seeing that the anova is approach is slightly... $\endgroup$ Commented Jul 9, 2018 at 15:58
-
$\begingroup$ ... underpowered relative to regression (p-values from anova are systematically higher). 3) The simulation you used a total sample size of 50, or on average, say, 17 observations per group. From what I'm seeing, the p-values from the two approaches diverge more strongly as the number of observations per group becomes smaller, becoming notable maybe when n-per-group is less than 10. 4) Why not just use logistic regression? $\endgroup$ Commented Jul 9, 2018 at 16:08
-
$\begingroup$ Thanks for the comment. I modified the comment in the code to reflect your correction about power. On the second comment, I only set out to show that it behaves as expected under the null. I typed and ran the code on a phone, so it wasn't as extensive as I would have liked to. In my simulations, the power of both methods track very closely. $\endgroup$ Commented Jul 9, 2018 at 16:09
-
$\begingroup$ At such small sample sizes, I wouldn't trust any estimates I obtained anyway especially in social science research where there is usually a lot of sample variance. Additionally, most effects are small such that for a well-powered study, the differences would probably not matter. My final reason would be the simplicity of one model relative to the other. $\endgroup$ Commented Jul 9, 2018 at 16:13
-
$\begingroup$ Out of curiosity, it would be interesting to know if a weighted least squares approach would match the power of logistic regression at small n. Weighting by the inverse of the standard deviation. In samples, $\sqrt{\hat{y}(1-\hat{y})}$. $\endgroup$ Commented Jul 9, 2018 at 16:21
$\begingroup$
$\endgroup$
No, ANOVA is not the appropriate test. You will want to conduct a chi-square test-of-independence test. This will involve generating a table that has counts for the 6 scenarios: treatment (3 options) by outcome (2 options).
