From: Austin, Peter <peter.austin_at_ices.on.ca>

Date: Tue 28 Feb 2006 - 05:46:26 EST

# Use a random effects logistic regression model with cluster-specific

# random intercepts.

resp.mat <- cbind(y.j,x.j)

# First column is number of successes, second column is number of failures.

# One row per cluster.

# Indicator variable for which group the cluster belongs to.

This email may contain confidential and/or privileged inform...{{dropped}}

R-help@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 Received on Tue Feb 28 05:56:58 2006

Date: Tue 28 Feb 2006 - 05:46:26 EST

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@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 Received on Tue Feb 28 05:56:58 2006

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:42:46 EST
*