Simulate data for stratified Cox model

I need to generate simulated data for a conditional logistic regression using R.

The code I wrote appears correct, but something is wrong. After generating my data, betas, and labels, I use glm() to attempt to estimate the betas. If everything is correct, the estimated betas from glm should be very close to the true betas. However, glm() fails to converge.

Below is my code. Can anyone point out what I'm doing wrong?

Thanks!

library(survival)

# Initialize parameters
k <- 5
n <- 10
g <- 20

p <- rep(0, n*g)
y <- rep(0, n*g)

# Create the group index
GI <- vector(mode="numeric")
for(group in 1:g){
GI <- c(GI, rep(group,n))
}

# generate the data
betas <- matrix(rnorm(k), ncol=1)
x <- matrix(rnorm(n*g*k), ncol=k)

# calculate the true labels
for(group in 1:g){
p[GI==group] <- exp(x[GI==group,] %*% betas) / sum( exp(x[GI==group,] %*% betas) )
y[n*(group-1)+which.max(p[GI==group])] <- 1
}

cl <- clogit(y ~x + strata(GI))

-
If you increase the variance in x or the variance in betas it converges. I usually also set iter.max to 10000 or so. –  Seth Apr 24 '12 at 6:03
Seth. Why would increasing variance affect convergence? –  Noah Apr 24 '12 at 8:57
I honestly don't know the particulars of why it works. Wish I had a better answer. –  Seth Apr 24 '12 at 14:27
@Noah As Seth asked below, why are you drawing random betas as well? –  Jason Morgan Jun 26 '12 at 1:56

I am not sure if you need your variance of x and beta to be 1? If you set the variance a bit higher your syntax will converge. change either:

betas <- matrix(rnorm(k,0,3), ncol=1)


or

x <- matrix(rnorm(n*g*k,0,3), ncol=k)


This trick works in many circumstances when the optimization routine is failing.

If you would like to know why, you could try plotting a couple profiles of the likelihood. Also, I am not sure why you are setting the betas to be random variables centered at zero. If you can adjust the means away from zero you may be able to converge with the variances closer to 1.

Good luck

-