[R] How to make try to catch warnings in logistic glm
Fredrik Nilsson
laf.nilsson at gmail.com
Mon Jun 22 14:55:32 CEST 2009
Dear list,
>From an earlier post I got the impression that one could promote
warnings from a glm to errors (presumably by putting
options(warn=1)?), then try() would flag them as errors. I’ve spent
half the day trying to do this, but no luck. Do you have an explicit
solution?
My problems is that I am trying to figure out during what conditions
one may find 5 significant parameters in a logistic regression with
only 21 bad and 28 outcomes (this appears in a published article) and
I would like to have my.glm<-try(glm(outcome~ald + fat+ hstop + btun
+ time, data=DF, family=binomial)) to indicate if warnings
(convergence / probabilities equal to 0/1) occurs (see attempt to code
below if my explanation is terse).
Best regards,
Fredrik Nilsson
------------------------------------------------------------------------------
Fredrik Nilsson, PhD
Competence Centre for Clinical Research
Lund University Hospital
Warnings I'd like to catch:
Warning in glm.fit(x = X, y = Y, weights = weights, start = start,
etastart = etastart, :
algorithm did not converge
Warning in glm.fit(x = X, y = Y, weights = weights, start = start,
etastart = etastart, :
fitted probabilities numerically 0 or 1 occurred
script
ngood<-28
nbad<-21
resmat<-"DFr"
j<-0
nsim<-100
res<-matrix(NA, nrow=nsim, ncol=6)
ow <- options("warn")
options(warn= 1)
outcome<-rep(c("G","B"), times=c(ngood, nbad))
outcome<-as.factor(outcome)
agood<-rep(c(0,1), times=c(27,1))
abad<-rep(c(0,1), times=c(14,7))
thgood<-rep(c(0,1), times=c(24,4))
thbad<-rep(c(0,1), times=c(12,9))
for (i in 1:nsim)
{
tgood<-sample(agood, ngood)
tbad<-sample(abad, nbad)
hstop<-c(tgood, tbad)
hstop<-as.factor(hstop)
aldgood<-rnorm(ngood, mean=54, sd=8/0.675)
aldbad<-rnorm(nbad, mean=64, sd=8/0.675)
ald<-c(aldgood, aldbad)
fatgood<-exp(rnorm(ngood, mean=log(23.9)-0.06^2/2,sd=0.06))
fatbad<-exp(rnorm(nbad, mean=log(27.4)-0.09^2/2, sd=0.09))
fat<-c(fatgood, fatbad)
thgood<-sample(thgood, ngood)
thbad<-sample(thbad, nbad)
btun<-c(thgood, thbad)
btun<-as.factor(btun)
timegood<-exp(rnorm(ngood, mean=log(443)-0.5^2/2, sd=0.5))
timebad<- exp(rnorm(nbad, mean=log(555)-0.6^2/2, sd=0.6))
time<-c(timegood, timebad)
DF<-data.frame(outcome, ald, fat, hstop, btun, time)
my.glm<-try(glm(outcome~ald + fat+ hstop + btun + time, data=DF,
family=binomial))
test<-!is.null(attr(my.glm, "try-error"))
if(test | !my.glm$converged)
{
res[i, 1:6]<-NA
} else
{
W<-summary(my.glm)$coefficients
lres<-as.numeric(W[,4]<=0.05)
if (sum(lres)==5)
{
j<-j+1
assign(paste(resmat, j, sep=".", collapse=""), DF)
}
res[i,1:6]<-lres
}
}
options(ow)
From: Dieter Menne <dieter.menne_at_menne-biomed.de>
Date: Fri, 14 Dec 2007 20:49:20 +0000 (UTC) Caroline Paulsen <cpaulsen
<at> u.washington.edu> writes:
>
> I'm attempting to run 250 permutations of a negative binomial GLM
> model for data on fish counts. Many of the models are fit
> appropriately, but some issue warnings such as "convergence not
> reached" or "step size truncated due to divergence."
....
> I'd like to figure out a way to flag the models that have warnings and
> output them into either a table or into R console. Then I could
> evaluate the problems associated with these models and know not choose
> them as the best models for fitting the data.
You could promote the options to errors (?warning) and use try() Dieter
More information about the R-help
mailing list