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 at csiro.au
http://www.cmis.csiro.au/bill.venables/
-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Roland R Regoes
Sent: Sunday, 15 January 2006 11:04 PM
To: r-help at 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 at 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