I would like to compare two logit models. One is a mixture of chemicals 'x' and 'y', the other is a model of just 'x'. I need help interpreting the 3 types of anova presented in the code/graphic below
However..
*A similar, but more complex situation: I would eventually like to compare independent tests of chemicals x and y, with a mixture of x and y. objective: look for possible synergistic effects (whether mixture is more potent, than either isolated chemical)
Would averaging the independent tests (x+y)/2
, making a model from averaged values, and then comparing this to a model made from the results of a true mixture of xy be folly?
Additional /supporting information can be found at bottom of post
#I prefer to use drc to create the glm
library(drc)
mod2 <- drm(probability ~ (dose), weights = total, data = mydata2, type ="binomial", fct=LL.2(names=c("Slope:b","ED50:e")))
mod10 <- drm(probability ~ (dose), weights = total, data = mydata10, type ="binomial", fct=LL.2(names=c("Slope:b","ED50:e")))
anova(mod2,mod10) #uses anova built into drc package
anova(mod2, mod10, test="F")
anova(mod2, mod10, test="Chisq")
#calculating Odd's Ratios (not included in output graphic below)
exp(coef(mod2))
exp(coef(mod10))
output:
Data:
mydata10 <-structure(list(dose = c(0, 0, 0, 3, 3, 3, 10, 10, 10, 27.5,27.5, 27.5, 50, 50, 50, 82.5, 82.5, 82.5), total = c(25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L,25L, 25L), affected = c(1, 0, 2.7, 3, 5.2, 5.3, 6.6, 7.3, 10.4, 14.3, 10.1, 14.7, 15.9, 12, 19.8, 20.4, 15.1, 19.4), probability = c(0.04,
0, 0.108, 0.12, 0.208, 0.212, 0.264, 0.292, 0.416, 0.572, 0.404,
0.588, 0.636, 0.48, 0.792, 0.816, 0.604, 0.776)), .Names = c("dose",
"total", "affected", "probability"), row.names = c(NA, -18L), class = "data.frame")
mydata2 <- structure(list(dose = c(0L, 0L, 0L, 3L, 3L, 3L, 10L, 10L, 10L,
30L, 30L, 30L, 55L, 55L, 55L, 85L, 85L, 85L), total = c(25L,
25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L,
25L, 25L, 25L, 25L), affected = c(0, 0, 5.4, 4.2, 8.6, 9.6, 7.2,
13, 14, 17.2, 17.2, 16.8, 19.2, 15, 20.2, 21.8, 22.6, 16.2),
probability = c(0, 0, 0.216, 0.168, 0.344, 0.384, 0.288,
0.52, 0.56, 0.688, 0.688, 0.672, 0.768, 0.6, 0.808, 0.872,
0.904, 0.648)), .Names = c("dose", "total", "affected", "probability"
), row.names = c(NA, -18L), class = "data.frame")
Supporting info:
when comparing the models with anova.drc the log likelihoods are negative values, with a p of 0. When calculated with the argument anova(mod2, mod10, test="F")
, the RSS (?) are positive values and the p is .1322. Which test should I use? Or would it be better to just compare the odd's ratios or some other method, e.g. (exp(coef(mod2)) - exp(coef(mod10))
?
Plos1 article from drc author: Supporting Info Describes the use drc's anova with some examples
From the anova.drc Cran page it seems anova.drc should only be used to test nested models, but I am not certain of this - in the results below you can see it produces the same result as anova(mod2, mod10, test="Chisq").
From anova.drc Cran page"
Specifying only a single object gives a test for lack-of-fit, comparing the non-linear regression model to a more general one-way or two-way ANOVA model. If two objects are specified a test for reduction from the larger to the smaller model is given. (This only makes statistical sense if the models are nested, that is: one model is a submodel of the other model.)
stackoverflow recommends against use of r^2
stats.stackexchange thread recommended anova,test="Chisq" to compare 2 glm models
Apologies for strange formatting, Unable to post more than 2 links due to no reputation, apologies - I normally just use stackoverflow.