1
$\begingroup$

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:

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.

$\endgroup$

1 Answer 1

1
$\begingroup$

"I need help interpreting the 3 types of anova presented in the code/graphic below"

The test that you are trying to do is generally used for model selection on the same data set (i.e either mod2 or mod10). A significant p-value means you should not simplify your model from a 4 -parameter logistic model (LL.4) to a 2-parameter logistic model (LL.2). Because you've set the data as "Binomial", I'd imagine the default Chi-squared test would be appropriate.

" 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)"

This can be done in various ways. If you had a fully-crossed design, you could fit a GLM and look for an interaction. However, in your example you don't. A simple way to see if there is a difference is to compare ED50 values and see if the 95% CI overlap.

mod2.ed50=ED(mod2, 50, interval = 'delta')
mod10.ed50=ED(mod10, 50, interval = 'delta')

mod10.ed50 is sig. greater than mod2.ed50. Also look into Concentration Addition models.

A final note. You refer to you models as GLMs but I am not sure this is the case as I get a different model fit using the regular glm function. It might be that adding 'type="binomial" to the drm model just allows weighting of the data. I stand to be corrected.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.