[R] Finding LD50 from an interaction Generalised Linear model
Bill.Venables at csiro.au
Bill.Venables at csiro.au
Wed Feb 13 01:36:19 CET 2008
The trick is to fit the model in a form which has the two separate
intercepts and the two separate slopes as the parameters.
You do have to realise that a*b, a*b-1, a/b, a/b-1, ... all specify the
same model, they just use different parameters for the job. (Yes,
really!)
> dat
ldose sex numdead
1 0 M 0
2 1 M 3
3 2 M 9
4 3 M 16
5 4 M 18
6 5 M 20
7 0 F 0
8 1 F 2
9 2 F 6
10 3 F 10
11 4 F 11
12 5 F 14
> dat <- transform(dat, Tot = 20)
> fm <- glm(numdead/20 ~ sex/ldose-1, binomial, dat, weights = Tot)
> coef(fm)
sexF sexM sexF:ldose sexM:ldose
-2.7634338 -3.4853625 0.7793144 1.5877754
> dose.p(fm, c(1,3)) ## females
Dose SE
p = 0.5: 3.545981 0.3025148
> dose.p(fm, c(2,4)) ## males
Dose SE
p = 0.5: 2.195123 0.1790317
>
In fact if you look at the book from which dose.p comes, you will find
an example not unlike this one done this way. At least I think that's
what I, er, read.
W.
Bill Venables
CSIRO Laboratories
PO Box 120, Cleveland, 4163
AUSTRALIA
Office Phone (email preferred): +61 7 3826 7251
Fax (if absolutely necessary): +61 7 3826 7304
Mobile: +61 4 8819 4402
Home Phone: +61 7 3286 7700
mailto:Bill.Venables at csiro.au
http://www.cmis.csiro.au/bill.venables/
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of Greme
Sent: Wednesday, 13 February 2008 3:17 AM
To: r-help at r-project.org
Subject: [R] Finding LD50 from an interaction Generalised Linear model
Hi,
I have recently been attempting to find the LD50 from two predicted fits
(For male and females) in a Generalised linear model which models the
effect
of both sex + logdose (and sex*logdose interaction) on proportion
survival
(formula = y ~ ldose * sex, family = "binomial", data = dat (y is the
survival data)). I can obtain the LD50 for females using the dose.p()
command in the MASS library with dose.p(mod1,c(1,2)). However I cannot
find
a way to determine the LD50 of males.
Any help on finding this male LD50 would be appreciated.
Pasting of R workspace below:
> rm(list=ls())
>
> ##checking file
> dat
ldose sex numdead
1 0 M 0
2 1 M 3
3 2 M 9
4 3 M 16
5 4 M 18
6 5 M 20
7 0 F 0
8 1 F 2
9 2 F 6
10 3 F 10
11 4 F 11
12 5 F 14
> str(dat)
'data.frame': 12 obs. of 3 variables:
$ ldose : int 0 1 2 3 4 5 0 1 2 3 ...
$ sex : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 1 1 1 1 ...
$ numdead: int 0 3 9 16 18 20 0 2 6 10 ...
> ##convert numdead to propdead
> dat$propdead<-dat$numdead/20
> ##Calculate survival from dead
> dat$numsurv<-20-dat$numdead
> dat$propsurv<-dat$numsurv/20
> ##check table
> dat
ldose sex numdead propdead numsurv propsurv
1 0 M 0 0.00 20 1.00
2 1 M 3 0.15 17 0.85
3 2 M 9 0.45 11 0.55
4 3 M 16 0.80 4 0.20
5 4 M 18 0.90 2 0.10
6 5 M 20 1.00 0 0.00
7 0 F 0 0.00 20 1.00
8 1 F 2 0.10 18 0.90
9 2 F 6 0.30 14 0.70
10 3 F 10 0.50 10 0.50
11 4 F 11 0.55 9 0.45
12 5 F 14 0.70 6 0.30
> str(dat)
'data.frame': 12 obs. of 6 variables:
$ ldose : int 0 1 2 3 4 5 0 1 2 3 ...
$ sex : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 1 1 1 1 ...
$ numdead : int 0 3 9 16 18 20 0 2 6 10 ...
$ propdead: num 0 0.15 0.45 0.8 0.9 1 0 0.1 0.3 0.5 ...
$ numsurv : num 20 17 11 4 2 0 20 18 14 10 ...
$ propsurv: num 1 0.85 0.55 0.2 0.1 0 1 0.9 0.7 0.5 ...
> ##load the lattice library
> library(lattice)
> ##plot data with mag + line width 1.5
> xyplot(propsurv~ldose,groups=sex,data=dat,cex=1.5,lwd=1.5,type="b")
> ##create model
> dat$n<-c(20)
> y<-cbind(dat$numsurv,dat$n-dat$numsurv)
> y
[,1] [,2]
[1,] 20 0
[2,] 17 3
[3,] 11 9
[4,] 4 16
[5,] 2 18
[6,] 0 20
[7,] 20 0
[8,] 18 2
[9,] 14 6
[10,] 10 10
[11,] 9 11
[12,] 6 14
> mod1<-glm(y~ldose*sex,dat,family="binomial")
> summary(mod1)
Call:
glm(formula = y ~ ldose * sex, family = "binomial", data = dat)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.94787 -0.36158 0.04914 0.63592 1.56417
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.7634 0.5231 5.282 1.28e-07 ***
ldose -0.7793 0.1550 -5.028 4.96e-07 ***
sexM 0.7219 0.8477 0.852 0.39444
ldose:sexM -0.8085 0.3131 -2.582 0.00981 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 136.1139 on 11 degrees of freedom
Residual deviance: 6.8938 on 8 degrees of freedom
AIC: 42.794
Number of Fisher Scoring iterations: 4
> anova(mod1,test="Chisq")
Analysis of Deviance Table
Model: binomial, link: logit
Response: y
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 11 136.114
ldose 1 106.323 10 29.791 6.266e-25
sex 1 15.063 9 14.729 1.040e-04
ldose:sex 1 7.835 8 6.894 0.005
> ##plot data
> pr<-expand.grid(sex=levels(dat$sex),ldose=seq(0,5,0.2))
> pr2<-data.frame(pr,preds=predict(mod1,type="response",newdata=pr))
>
> mm<-dat[dat$sex=="M",]
> ff<-dat[dat$sex=="F",]
>
> par(mfrow=c(1,1))
>
plot(mm$numsurv/mm$n~ldose,mm,cex=1.5,cex.axis=1.5,xlab="logDose",ylab="
Proportion
> Survived")
> points(ff$numsurv/ff$n~ldose,ff,cex=1.5,col=2)
>
> ##prediction lines
> lines(pr2[pr2$sex=="M",]$preds~pr2[pr2$sex=="M",]$ldose)
> lines(pr2[pr2$sex=="F",]$preds~pr2[pr2$sex=="F",]$ldose,col=2)
>
> ##lethaldose50line+legend
> abline(h=0.5,lty=2)
> text(0.5,0.48,"Find LD50")
> legend(3.5,1,c("Male","Female"),col=1:2,lty=1,cex=1.5)
> library(MASS)
> dose.p(mod1,c(1,2))
Dose SE
p = 0.5: 3.545981 0.3025148
>
--
View this message in context:
http://www.nabble.com/Finding-LD50-from-an-interaction-Generalised-Linea
r-model-tp15436597p15436597.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list