R help - Jun 2009 - How to make try to catch warnings in logistic glm

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

Fredrik Nilsson-5 wrote:
> 
>>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?
> 
> 
You were rather close... 

1) use warn=2
2) make sure to get the | and || right (this is nasty, I never learn it)

Dieter


 options(warn=2) #
 my.glm<-try(glm(outcome~ald + fat+ hstop + btun + time, data=DF,
family=binomial),FALSE)
 failed <- inherits(my.glm,"try-error")
 if(failed || !my.glm$converged)
 {
   res[i, 1:6]<-NA
 } else
 {


-- 
View this message in context:
http://www.nabble.com/How-to-make-try-to-catch-warnings-in-logistic-glm-tp24146726p24147228.html
Sent from the R help mailing list archive at Nabble.com.