I am trying to do a dose response in my dataset, but nothing go a head.
I am adapting a script shared on the web, but I unable to make it
useful for my dataset. I would like to got the LC50 for each Isolado
and if there are differences between then.
My data is https://dl.dropboxusercontent.com/u/34009642/R/dead_alive.csv
Here what I copy and try to modifying:
library(plyr)
library(lattice)
library(lme4)
library(arm)
library(lmerTest)
library(faraway)
library(car)
## Conc are concentration. I input only the coef, but, all,
## except 0, that is my control (without Isolado), are base
## 10. i.e: 10^4, 10^6 e 10^8.
data <- read.table("dead_alive.csv", sep = "\t",
dec=",", header = TRUE)
data$Rep <- factor(data$Rep)
mean_data <- ddply(data, c("Isolado", "Conc",
"Day"), numcolwise(mean))
xyplot(Dead/(Dead + Live) ~ Conc|Isolado, groups = Day, type = "l",
ylab='Probability', xlab='Dose', data = mean_data)
xyplot(Dead/(Dead + Live) ~ Day|Isolado, groups = Conc, type = "l",
ylab='Probability', xlab='Dose', data = mean_data)
model.logit <- glmer(cbind(Dead, Live) ~ -1 + Isolado + Isolado:Conc +
(0 + Conc|Day), family=binomial, data = data)
Anova(model.logit)
summary(model.logit)
model.probit <- glmer(cbind(Dead, Live) ~ Isolado + Isolado:Conc + (0
+ Conc|Day), family=binomial(link=probit), data=data)
model.cloglog <- glm(cbind(Dead, Live) ~ Isolado + Isolado:Conc + (1 +
Conc|Day), family=binomial(link=cloglog), data=data)
x <- seq(0,8, by=0.2)
prob.logit <- ilogit(model.logit$coef[1] + model.logit$coef[2]*x)
prob.probit <- pnorm(model.probit$coef[1] + model.probit$coef[2]*x)
prob.cloglog <- 1-exp(-exp((model.cloglog$coef[1] +
model.cloglog$coef[2]*x)))
with(subdata, plot(Dead/(Dead + Live) ~ Conc, group = Day, )
lines(x, prob.logit) # solid curve = logit
lines(x, prob.probit, lty=2) # dashed = probit
lines(x, prob.cloglog, lty=5) # longdash = c-log-log
plot(x, prob.logit, type='l', ylab='Probability',
xlab='Dose') # solid
curve = logit
lines(x, prob.probit, lty=2) # dashed = probit
lines(x, prob.cloglog, lty=5) # longdash = c-log-log
matplot(x, cbind(prob.probit/prob.logit,
(1-prob.probit)/(1-prob.logit)), type='l', xlab='Dose',
ylab='Ratio')
matplot(x, cbind(prob.cloglog/prob.logit,
(1-prob.cloglog)/(1-prob.logit)), type='l', xlab='Dose',
ylab='Ratio')
model.logit.data <- glm(cbind(Dead,Live) ~ Conc, family=binomial,
data=data)
pred2.5 <- predict(model.logit.data, newdata=data.frame(Conc=2.5), se=T)
ilogit(pred2.5$fit)
ilogit(c(pred2.5$fit - 1.96*pred2.5$se.fit, pred2.5$fit +
1.96*pred2.5$se.fit))
## what are this 1.96???? Where it come from?
### If there are several predictors, just put in the code
### above something like:
### newdata=data.frame(conc=2.5,x2=4.6,x3=5.8)
### or whatever is the desired set of predictor values...
### Effective Dose calculation:
# What is the concentration that yields a probability of 0.5 of an
# insect dying?
library(MASS)
dose.p(model.logit.data, p=0.5)
# A 95% CI for the ED50:
c(2 - 1.96*0.1466921, 2 + 1.96*0.1466921)
# What is the concentration that yields a probability of 0.8 of an
# insect dying?
dose.p(model.logit.data, p=0.8)
--
Laia, ML