# Re: [R] Finding LD50 from an interaction Generalised Linear model

From: <Bill.Venables_at_csiro.au>
Date: Wed, 13 Feb 2008 10:36:19 +1000

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

```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
>

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

```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 ...
```

> dat\$propsurv<-dat\$numsurv/20
> ##check table
> dat

```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 ...
```

> 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

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
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
and provide commented, minimal, self-contained, reproducible code.
```
Received on Wed 13 Feb 2008 - 00:40:05 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 13 Feb 2008 - 10:30:13 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.