R help - Apr 2012 - Calculating NOEL using R and logistic regression - Toxicology

Hello, I used the glm function in R to fit a dose-response relationship and
then have been using dose.p to calculate the LC50, however I would like to
calculate the NOEL (no observed effect level), ie the lowest dose above
which responses start occurring. Does anyone know how to do this?

	[[alternative HTML version deleted]]

dear Danielle,

The NOEL is a threshold value or breakpoint in the range of dose. Have a look to
the
package segmented to estimate a GLM with unknown breakpoints. The code
(untested) should
be something like

library(segmented)
o<-glm(y~1, family=binomial)
os<-segmented(o, ~dose, psi=starting_psi)

Also the package segmented includes the dataset down that can be useful as an
example..

data(down)
with(down, plot(age, cases/births))

There is a paper of mine on R news 2008 discussing the package..

hope this helps you,
vito



On Mon, 2 Apr 2012 14:45:06 -0800, Danielle Duncan wrote
> Hello, I used the glm function in R to fit a dose-response relationship and
> then have been using dose.p to calculate the LC50, however I would like to
> calculate the NOEL (no observed effect level), ie the lowest dose above
> which responses start occurring. Does anyone know how to do this?
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> 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.

--
Open WebMail Project (http://openwebmail.org)

Is this the smallest observed dose that has an effect? If so, then you
don't need the glm to find it. Here is a simulated example:

set.seed(101)
X = rep(1:10,each=10)
lp = -5 + 0.5*X
Y = rbinom(length(X),size=1,p=1/(1+exp(-lp)))
# is this the NOEL?
min(X[Y==1]) # picks out observations with adverse effects, chooses the
smallest value
glmfit = glm(Y~X,family=binomial)
plot(1:10, predict(glmfit,newdata=data.frame(X=1:10),type="response"),
type="l",ylim=c(0,1),xlab="X",ylab="Y")
rug(jitter(X[Y==0]),side=1)
rug(jitter(X[Y==1]),side=3)

On Tue, Apr 3, 2012 at 3:19 PM, Danielle Duncan
<dlduncan2@alaska.edu>wrote:
> Thanks, that is interesting, but what I'm really after is an easy
"no
> observed effect level", using a binomial logistic model ie glm. Have a
> great day!
>
> On Mon, Apr 2, 2012 at 3:38 PM, vito.muggeo <vito.muggeo@unipa.it>
wrote:
>
> > dear Danielle,
> >
> > The NOEL is a threshold value or breakpoint in the range of dose. Have
a
> > look to the
> > package segmented to estimate a GLM with unknown breakpoints. The code
> > (untested) should
> > be something like
> >
> > library(segmented)
> > o<-glm(y~1, family=binomial)
> > os<-segmented(o, ~dose, psi=starting_psi)
> >
> > Also the package segmented includes the dataset down that can be
useful
> as
> > an example..
> >
> > data(down)
> > with(down, plot(age, cases/births))
> >
> > There is a paper of mine on R news 2008 discussing the package..
> >
> > hope this helps you,
> > vito
> >
> >
> >
> > On Mon, 2 Apr 2012 14:45:06 -0800, Danielle Duncan wrote
> > > Hello, I used the glm function in R to fit a dose-response
relationship
> > and
> > > then have been using dose.p to calculate the LC50, however I
would like
> > to
> > > calculate the NOEL (no observed effect level), ie the lowest dose
above
> > > which responses start occurring. Does anyone know how to do this?
> > >
> > >       [[alternative HTML version deleted]]
> > >
> > > ______________________________________________
> > > R-help@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.
> >
> >
> > --
> > Open WebMail Project (http://openwebmail.org)
> >
> >
>
>        [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help@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.
>


-- 
Drew Tyre

School of Natural Resources
University of Nebraska-Lincoln
416 Hardin Hall, East Campus
3310 Holdrege Street
Lincoln, NE 68583-0974

phone: +1 402 472 4054
fax: +1 402 472 2946
email: atyre2@unl.edu
http://snr.unl.edu/tyre
http://aminpractice.blogspot.com
http://www.flickr.com/photos/atiretoo

	[[alternative HTML version deleted]]

Thanks for the response, I should have clarified that the NOEL is the
smallest dose above which there is a statistically significant effect.

On Tue, Apr 3, 2012 at 12:44 PM, Drew Tyre <atyre2@unl.edu> wrote:
> Is this the smallest observed dose that has an effect? If so, then you
> don't need the glm to find it. Here is a simulated example:
>
> set.seed(101)
> X = rep(1:10,each=10)
> lp = -5 + 0.5*X
> Y = rbinom(length(X),size=1,p=1/(1+exp(-lp)))
> # is this the NOEL?
> min(X[Y==1]) # picks out observations with adverse effects, chooses the
> smallest value
> glmfit = glm(Y~X,family=binomial)
> plot(1:10,
predict(glmfit,newdata=data.frame(X=1:10),type="response"),
> type="l",ylim=c(0,1),xlab="X",ylab="Y")
> rug(jitter(X[Y==0]),side=1)
> rug(jitter(X[Y==1]),side=3)
>
> On Tue, Apr 3, 2012 at 3:19 PM, Danielle Duncan
<dlduncan2@alaska.edu>wrote:
>
>> Thanks, that is interesting, but what I'm really after is an easy
"no
>> observed effect level", using a binomial logistic model ie glm.
Have a
>> great day!
>>
>> On Mon, Apr 2, 2012 at 3:38 PM, vito.muggeo
<vito.muggeo@unipa.it> wrote:
>>
>> > dear Danielle,
>> >
>> > The NOEL is a threshold value or breakpoint in the range of dose.
Have a
>> > look to the
>> > package segmented to estimate a GLM with unknown breakpoints. The
code
>> > (untested) should
>> > be something like
>> >
>> > library(segmented)
>> > o<-glm(y~1, family=binomial)
>> > os<-segmented(o, ~dose, psi=starting_psi)
>> >
>> > Also the package segmented includes the dataset down that can be
useful
>> as
>> > an example..
>> >
>> > data(down)
>> > with(down, plot(age, cases/births))
>> >
>> > There is a paper of mine on R news 2008 discussing the package..
>> >
>> > hope this helps you,
>> > vito
>> >
>> >
>> >
>> > On Mon, 2 Apr 2012 14:45:06 -0800, Danielle Duncan wrote
>> > > Hello, I used the glm function in R to fit a dose-response
>> relationship
>> > and
>> > > then have been using dose.p to calculate the LC50, however I
would
>> like
>> > to
>> > > calculate the NOEL (no observed effect level), ie the lowest
dose
>> above
>> > > which responses start occurring. Does anyone know how to do
this?
>> > >
>> > >       [[alternative HTML version deleted]]
>> > >
>> > > ______________________________________________
>> > > R-help@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.
>> >
>> >
>> > --
>> > Open WebMail Project (http://openwebmail.org)
>> >
>> >
>>
>>        [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help@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.
>>
>
>
>
> --
> Drew Tyre
>
> School of Natural Resources
> University of Nebraska-Lincoln
> 416 Hardin Hall, East Campus
> 3310 Holdrege Street
> Lincoln, NE 68583-0974
>
> phone: +1 402 472 4054
> fax: +1 402 472 2946
> email: atyre2@unl.edu
> http://snr.unl.edu/tyre
> http://aminpractice.blogspot.com
> http://www.flickr.com/photos/atiretoo
>
	[[alternative HTML version deleted]]

Dear Danielle,

At least in industrial toxicology (my original background) the recent
tendency has been to use benchmark dose (BSD) approach instead of NOEL
or NOAEL approach due to various problems with the definition and
estimation of NO(A)EL. In R this can be achieved using the packages
drc and bmd as already mentioned by Drew Tyre.

Here's a short code example that gives you the BMD at 1% level:

library(drc)
library(bmd)
data(daphnids)
d24<-daphnids[daphnids$time=="24h",]
fit <- drm(no/total~dose, weights = total, data = d24, fct = LL.2(),
type = "binomial")
ED(fit, 1)
bmd(FIT, 0.01)

See the following documents for more information, if you're interested
in using BSD instead of NO(A)EL:

http://www.cdpr.ca.gov/docs/risk/bmdquant.pdf
http://www.efsa.europa.eu/en/efsajournal/doc/1150.pdf

Jarno

Make that

bmd(fit, 0.01)

in my previous post.

Jarno

Thanks everyone for the advice, you raise interesting points. Maybe the
best thing for me to do is do an ANOVA in R with binomial data (if
possible) and find the lowest dose that gives a significant difference from
the controls.

On Mon, Apr 2, 2012 at 2:45 PM, Danielle Duncan
<dlduncan2@alaska.edu>wrote:
> Hello, I used the glm function in R to fit a dose-response relationship
> and then have been using dose.p to calculate the LC50, however I would like
> to calculate the NOEL (no observed effect level), ie the lowest dose above
> which responses start occurring. Does anyone know how to do this?



-- 
"Wilderness isn’t the wide open spaces. It’s the wild things in it”

	[[alternative HTML version deleted]]

I suppose I'll just report a LC10 using the dose.p function in the package
MASS using my glm fitted logistic regression on binomial data. Thanks
everyone for ideas & input! The LOEC seems to be a flawed
calculation.......I'll research it. Again, thanks all!

On Mon, Apr 2, 2012 at 2:45 PM, Danielle Duncan
<dlduncan2@alaska.edu>wrote:
> Hello, I used the glm function in R to fit a dose-response relationship
> and then have been using dose.p to calculate the LC50, however I would like
> to calculate the NOEL (no observed effect level), ie the lowest dose above
> which responses start occurring. Does anyone know how to do this?



-- 
"Wilderness isn’t the wide open spaces. It’s the wild things in it”

	[[alternative HTML version deleted]]