Re: [R] problems with glm

From: Roland R Regoes <roland.regoes_at_env.ethz.ch>
Date: Mon 16 Jan 2006 - 00:56:13 EST

        I am sorry, I should have emphasized that: I really mean 'family = binomial(link="log")'.

	I am fitting infection data. Failure and success correspond to being infected or not. The probability of success (ie, not being infected in my context) is derived from a population model describing the interaction between hosts and parasites: 
	P(not infected) = exp(-InfectionRate*ParasiteDose)
'ParasiteDose' corresponds to 'predictor'. I want to estimate the infection rate. Hence I need 'binomial(link="log")' and must set intercept=0.

Roland

On Mon, Jan 16, 2006 at 12:16:09AM +1100, Bill.Venables@csiro.au wrote:
> Most of your problems seem to come from 'link = "log"' whereas you
> probably mean 'link = logit' (which is the default. Hence:
>
> ##########################################
> > success <- c(13,12,11,14,14,11,13,11,12)
> > failure <- c(0,0,0,0,0,0,0,2,2)
> > predictor <- c(0,80*5^(0:7))
> > glm(cbind(success,failure) ~ predictor,
> + family = binomial, #(link="log"),
> + control = glm.control(epsilon=1e-8,trace=TRUE,maxit=50))
> Deviance = 7.621991 Iterations - 1
> Deviance = 6.970934 Iterations - 2
> Deviance = 6.941054 Iterations - 3
> Deviance = 6.940945 Iterations - 4
> Deviance = 6.940945 Iterations - 5
>
> Call: glm(formula = cbind(success, failure) ~ predictor, family =
> binomial, control = glm.control(epsilon = 1e-08, trace = TRUE,
> maxit = 50))
>
> Coefficients:
> (Intercept) predictor
> 4.180e+00 -4.106e-07
>
> Degrees of Freedom: 8 Total (i.e. Null); 7 Residual
> Null Deviance: 12.08
> Residual Deviance: 6.941 AIC: 15.85
> >
>
> #######################################
>
> Bill Venables,
> CMIS, CSIRO Laboratories,
> PO Box 120, Cleveland, Qld. 4163
> AUSTRALIA
> Office Phone (email preferred): +61 7 3826 7251
> Fax (if absolutely necessary): +61 7 3826 7304
> Mobile (rarely used): +61 4 1963 4642
> Home Phone: +61 7 3286 7700
> mailto:Bill.Venables@csiro.au
> http://www.cmis.csiro.au/bill.venables/
>
>
>
> -----Original Message-----
> From: r-help-bounces@stat.math.ethz.ch
> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Roland R Regoes
> Sent: Sunday, 15 January 2006 11:04 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] problems with glm
>
>
> Dear R users,
>
> I am having some problems with glm. The first is an error message
> "subscript out of bounds". The second is the fact that reasonable
> starting values are not accepted by the function.
>
> To be more specific, here is an example:
>
> > success <- c(13,12,11,14,14,11,13,11,12)
> > failure <- c(0,0,0,0,0,0,0,2,2)
> > predictor <- c(0,80*5^(0:7))
> > glm(cbind(success,failure) ~ predictor + 0,
> + family=binomial(link="log"),
> + control=glm.control(epsilon=1e-8,trace=TRUE,maxit=50))
> Deviance = 3.348039 Iterations - 1
> Error: subscript out of bounds
> In addition: Warning message:
> step size truncated: out of bounds
> >
>
> The model with intercept yields:
>
> > glm(cbind(success,failure) ~ predictor ,
> + family=binomial(link="log"),
> + control=glm.control(epsilon=1e-8,trace=FALSE,maxit=50))
>
> Call: glm(formula = cbind(success, failure) ~ predictor, family =
> binomial(link = "log"), control = glm.control(epsilon = 1e-08,
> trace = FALSE, maxit = 50))
>
> Coefficients:
> (Intercept) predictor
> -5.830e-17 -4.000e-08
>
> Degrees of Freedom: 8 Total (i.e. Null); 7 Residual
> Null Deviance: 12.08
> Residual Deviance: 2.889 AIC: 11.8
> There were 33 warnings (use warnings() to see them)
> > warnings()
> 1: step size truncated: out of bounds
> ...
> 31: step size truncated: out of bounds
> 32: algorithm stopped at boundary value in: glm.fit(x = X, y = Y,
> weights = weights, start = start, etastart = etastart, ...
> 33: fitted probabilities numerically 0 or 1 occurred in: glm.fit(x = X,
> y = Y, weights = weights, start = start, etastart = etastart, ...
> >
>
> Since the intercept in the above fit is fairly small I thought I
> could use -4e-8 as a reasonable starting value in a model without
> intercept. But to no avail:
>
> > glm(cbind(success,failure) ~ predictor + 0, start=-4e-8,
> + family=binomial(link="log"),
> + control=glm.control(epsilon=1e-8,trace=FALSE,maxit=50))
> Error in glm.fit(x = X, y = Y, weights = weights, start = start,
> etastart = etastart, :
> cannot find valid starting values: please specify some
> >
>
> I am stuck here. Am I doing something wrong when specifying the
> starting value? I would appreciate any help. (I could not find anything
> relevant in the documentation of glm and the mailing list archives, but
> I did not read the source code of glm yet.)
>
> Roland
>
>
> PS: I am using R Version 2.2.0 (R Cocoa GUI 1.13 (1915)) on MacOSX
> 10.4.4
>
> --
> -----------------------------------------------------------------------
> Roland Regoes
> Theoretical Biology
> Universitaetsstr. 16
> ETH Zentrum, CHN H76.1
> CH-8092 Zurich, Switzerland
>
> tel: +41-44-632-6935
> fax: +41-44-632-1271
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html

-- 
-----------------------------------------------------------------------
Roland Regoes
Theoretical Biology
Universitaetsstr. 16 
ETH Zentrum, CHN H76.1 
CH-8092 Zurich, Switzerland

tel: +41-44-632-6935
fax: +41-44-632-1271

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Mon Jan 16 01:07:08 2006

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:42:04 EST