Folks,
There appears to be a bug(?) in step() when used to screen logistic
models. The problem appears to be specific to 1.2.1 (or maybe also
1.2.0?), as the proper behavior was observed in earlier versions. When
I did the same on the surrogate log linear model, everything seemed
okey.
The data involved was the detergent data found in earlier editions of
MASS, as given below,
> detg1
Temp M.user Soft M X
1 Low N Hard 42 68
3 High N Hard 30 42
5 Low Y Hard 52 37
7 High Y Hard 43 24
9 Low N Medium 50 66
11 High N Medium 23 33
13 Low Y Medium 55 47
15 High Y Medium 47 23
17 Low N Soft 53 63
19 High N Soft 27 29
21 Low Y Soft 49 57
23 High Y Soft 29 19
A constant model was fit to the data
> detg1.m0 <- glm(cbind(X,M)~1,binomial,detg1)
> detg1.m0
Call: glm(formula = cbind(X, M) ~ 1, family = binomial, data = detg1)
Coefficients:
(Intercept)
0.01587
Degrees of Freedom: 11 Total (i.e. Null); 11 Residual
Null Deviance: 32.83
Residual Deviance: 32.83 AIC: 92.52
The following was the results I got from R.1.2.1 on both AIX and FreeBSD,
> step(detg1.m0,scope=list(upper=~M.user*Temp*Soft))
Start: AIC= 92.52
cbind(X, M) ~ 1
Df Deviance AIC
<none> 32.83 92.52
+ M.user 1 1059.98 1121.67
+ Temp 1 2236.57 2298.27
+ Soft 2 2565.93 2629.63
Call: glm(formula = cbind(X, M) ~ 1, family = binomial, data = detg1)
Coefficients:
(Intercept)
0.01587
Degrees of Freedom: 11 Total (i.e. Null); 11 Residual
Null Deviance: 32.83
Residual Deviance: 32.83 AIC: 92.52
I then tried stepAIC() in the MASS package, although I was not sure
why it should be any different. The following was the results,
> stepAIC(detg1.m0,scope=list(lower=~.,upper=~.+M.user))
Start: AIC= 92.52
cbind(X, M) ~ 1
Df Deviance AIC
+ M.user 1 12.244 73.942
<none> 32.826 92.524
Step: AIC= 73.94
cbind(X, M) ~ M.user
Df Deviance AIC
- M.user 1 0.428 60.125
<none> 12.244 73.942
Step: AIC= 92.52
cbind(X, M) ~ 1
Call: glm(formula = cbind(X, M) ~ 1, family = binomial, data = detg1)
Coefficients:
(Intercept)
0.01587
Degrees of Freedom: 11 Total (i.e. Null); 11 Residual
Null Deviance: 32.83
Residual Deviance: 32.83 AIC: 92.52
Warning message:
non-integer #successes in a binomial glm! in: eval(expr, envir, enclos)
The AIC obtained by stepAIC() appeared to be correct in the first step
add1() operation, but the drop1() screwed up in the next step and
brought me right back to the constant model.
Will check on my linux box tonight, but this one might not be
platform-specific as I was able to duplicate on two platforms so far.
Chong
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Fri, 16 Feb 2001, Chong Gu wrote:> There appears to be a bug(?) in step() when used to screen logistic > models. The problem appears to be specific to 1.2.1 (or maybe also > 1.2.0?), as the proper behavior was observed in earlier versions. When > I did the same on the surrogate log linear model, everything seemed > okey.The bug is in add1.glm. It was intriduced in 1.2.0. It applies only to binomial models specified in cbind(success, failure) form. The problem is that what glm returns for the weights and y is in the 0 <= y <= 1 form, and that is not what model.extract gives. I will commit a fix shortly. stepAIC is wrong because of the same bug, but the factors come out different. -- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
perhaps this should be directed to r-devel or somesuch, but I thought
that others might wish to comment on the concept.
As a fairly new (and very enthusiastic and very appreciative) user of R I
find myself frequently going to the help files, quite often just for some
simple help on syntax or function arguments that I vaguely remember ..
hist(), for instance, has a number of options some of which are to do with
left and right boundaries of the cells that I sometimes need but not
often.
Some form of "code completion" might be of assistance to users like me
..
not exactly casual users or complete novices, but those who tend to use R
in bursts perhaps separated by weeks or longer. Code completion is
implemented in Delphi and perhaps in other languages that have a GUI IDE
(integrated development environment).. essentially it comes down to the
IDE prompting you with a function prototype after some delay or on
pressing some special key (eg ctrl+?)
Could something similar be done in R and is it desirable? .. there perhaps
might be some simple approximation to code completion achievable by some
modest changes to the syntax .. eg hist(?) might paste a full function
prototype like
hist(x, breaks, freq = NULL, probability = !freq,
include.lowest = TRUE,
right = TRUE, col = NULL, border = par("fg"),
main = paste("Histogram of" , xname),
xlim = range(breaks), ylim = NULL,
xlab = xname, ylab,
axes = TRUE, plot = TRUE, labels = FALSE,
nclass = NULL, ...)
into the console window (or perhaps some special edit window)
well, something to ponder, fwiw.
John Aitchison
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._