[R] singular convergence in glmmPQL

[Berton Gunter]

If I understand you correctly (your post was, ummm... a bit long),

?try

should help. Condition on the class of the return inheriting from
'try-error' . i.e.

result<-try(glimmPQL(...))
if(inherits(result,'try-error'))...

As the try man page says, fancier things can be done via tryCatch()

-- Bert Gunter
Genentech Non-Clinical Statistics
South San Francisco, CA
 
"The business of the statistician is to catalyze the scientific learning
process."  - George E. P. Box
 
 

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Austin, Peter
> Sent: Monday, February 27, 2006 10:46 AM
> To: r-help at stat.math.ethz.ch
> Cc: Austin, Peter
> Subject: [R] singular convergence in glmmPQL
> 
> 
> I am using the 'glmmPQL function in the 'MASS' library to fit 
> a mixed effects logistic regression model to simulated data.  
> I am conducting a series of simulations, and with certain 
> simulated datasets, estimation of the random effects logistic 
> regression model unexpectedly terminates.  I receive the 
> following error message from R:
> 
> Error in lme.formula(fixed=zz + arm.long,random=~1 | id.long, 
> method="ML", : singular convergence (7)
> 
> I would like to modify my code so that the resultant model is 
> assigned a null value, rather than having estimation terminate.
> 
> My R code for this example is contained below.
> 
> # Data are read in for the two sets of clusters.  There are 
> 100 clusters
> # in each of the two groups.
> 
> y.1j 
> <-c(3,2,3,1,3,3,2,2,3,1,1,3,2,2,3,3,3,1,5,2,4,4,1,3,3,4,5,3,3,
> 0,2,2,0,2,2,2,1,3,2,3,1,0,4,4,4,2,1,2,4,1,0,1,0,1,2,4,2,1,1,2,
> 3,3,1,3,0,2,3,4,2,2,4,2,1,3,1,4,3,0,3,4,3,0,1,3,1,0,2,2,6,2,2,
> 1,1,1,2,1,0,2,2,2)
> 
> y.2j <- 
> c(3,4,3,3,3,4,1,8,2,3,3,1,3,2,1,2,4,6,1,3,6,3,4,3,5,5,4,6,3,3,
> 4,0,4,2,2,4,1,0,5,2,7,4,4,3,3,4,4,6,1,1,2,2,2,4,0,2,1,5,5,2,3,
> 4,1,3,1,1,1,4,3,3,3,2,1,3,1,3,2,3,2,3,4,2,3,2,7,2,2,2,3,2,6,2,
> 2,3,3,3,2,3,1,4)
> # Number of positive responses in each cluster in each group.
> 
> n.clusters.1 <- length(y.1j)
> n.clusters.2 <- length(y.2j)
> # Number of clusters in each group.
> 
> m.1j <- rep(20,n.clusters.1)
> m.2j <- rep(20,n.clusters.2)
> # 20 subjects in each cluster in each group.
> 
> M.1 <- sum(m.1j)
> M.2 <- sum(m.2j)
> # Number of subjects in each of the two groups.
> 
> K <- n.clusters.1 + n.clusters.2
> # Total number of clusters in the two groups combined.
> 
> ##############################################################
> ##############
> # Use a random effects logistic regression model with cluster-specific
> # random intercepts. 
> ##############################################################
> ##############
> 
> library("MASS")
> # import MASS library for using function 'glmmPQL'
> 
> # create subject-specific responses from cluster-aggregate responses.
> y.j <- c(y.1j,y.2j)
> m.j <- c(m.1j,m.2j)
> x.j <- m.j - y.j
> resp.mat <- cbind(y.j,x.j)
> # First column is number of successes, second column is 
> number of failures.
> # One row per cluster.
> 
> y <- rep(c(1,0),K)
> pat.mat <- matrix(t(resp.mat),ncol=1,byrow=T)
> y.long <- rep(y,pat.mat)
> # Create vector of 1/0 outcomes with 1 observation per subject.
> 
> id.long <- rep(1:K,c(m.1j,m.2j))
> 
> arm.0 <- rep(1,M.1)
> arm.1 <- rep(0,M.2)
> arm.long <- c(arm.0,arm.1)
> # Indicator variable for which group the cluster belongs to.
> 
> re.1 <- glmmPQL(y.long ~ arm.long,random = ~1|id.long,family=binomial)
> 
> 
> My question is this:
> Is it possible for the vector "re.1" to be returned with a 
> null value if there is a singular convergence in the 
> underlying call to 'lme'? Rather than having the execution 
> terminate, can a null value be assigned to "re.1"?  As it 
> stands, "re.1" is undefined when estimation terminates.
> 
> Thanks very much,
> 
> Peter
> 
> Peter Austin, PhD
> Senior Scientist, Institute for Clinical Evaluative Sciences.
> Associate Professor, Departments of Public Health Sciences 
> and Health Policy, Management and Evaluation, University of Toronto.
>  
> Institute for Clinical Evaluative Sciences 
> G106 - 2075 Bayview Avenue 
> Toronto, Ontario, M4N 3M5 
> Tel:  (416) 480-6131
> Fax: (416) 480-6048 
> ICES Web Site: www.ices.on.ca
> 
> 
> ______________________________________________________________
> _____________ 
> This email may contain confidential and/or privileged 
> inform...{{dropped}}
> 
> ______________________________________________
> 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
>