[R] Dose response glmer

[Marcelo Laia]

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