[R] GLMM(..., family=binomial(link="cloglog"))?

[Spencer Graves]

	  Another GLMM/glmm problem:  I simulate rbinom(N, 100, pz), where 
logit(pz) = rnorm(N).  I'd like to estimate the mean and standard 
deviation of logit(pz).  I've tried GLMM{lme4}, glmmPQL{MASS}, and 
glmm{Jim Lindsey's repeated}.  In several replicates of this for N = 10, 
100, 500, etc., my glmm call produced estimates of the standard 
deviation of the random effect in the range of 0.6 to 0.8 (never as high 
as the 1 simulated).  Meanwhile, my calls to GLMM produced estimates 
between 1e-12 and 1e-9, while the glmmPQL results tended to be closer to 
0.001, though it gave one number as high as 0.7.  (I'm running R 1.9.1 
alpha, lme4 0.6-1 under Windows 2000)

	  Am I doing something wrong, or do these results suggest bugs in the 
software or deficiencies in the theory or ... ?

	  Consider the following:

 > set.seed(1); N <- 10
 > z <- rnorm(N)
 > pz <- inv.logit(z)
 > DF <- data.frame(z=z, pz=pz, y=rbinom(N, 100, pz)/100, n=100, 
smpl=factor(1:N))
 > GLMM(y~1, family=binomial, data=DF, random=~1|smpl, weights=n)
Generalized Linear Mixed Model

Fixed: y ~ 1
Data: DF
  log-likelihood:  -55.8861
Random effects:
  Groups Name        Variance   Std.Dev.
  smpl   (Intercept) 1.7500e-12 1.3229e-06

Estimated scale (compare to 1)  3.280753

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.148271   0.063419  2.3379  0.01939

Number of Observations: 10
Number of Groups: 10

 > library(repeated) # Jim Lindsey's repeated package
 > glmm(y~1, family=binomial, data=DF, weights=n, nest=smpl)

  Warning messages:
1: non-integer #successes in a binomial glm! in: eval(expr, envir, 
enclos) ...

Coefficients:
(Intercept)           sd
      0.1971       0.6909

Degrees of Freedom: 9 Total (i.e. Null);  8 Residual
Null Deviance:      111.8
Residual Deviance: 32.03        AIC: 1305
Normal mixing variance: 0.4773187

 > glmmPQL(y~1, family=binomial, data=DF, random=~1|smpl, weights=n)

Linear mixed-effects model fit by maximum likelihood
   Data: DF
   Log-likelihood: -10.78933
   Fixed: y ~ 1
(Intercept)
   0.1622299

Random effects:
  Formula: ~1 | smpl
         (Intercept)  Residual
StdDev:   0.7023349 0.5413203

Variance function:
  Structure: fixed weights
  Formula: ~invwt
Number of Observations: 10
Number of Groups: 10

	  Suggestions?

	  So far, the best tool I've got for this problem is a normal 
probability plot of a transform of the binomial responses with Monte 
Carlo confidence bands, as suggested by Venables and Ripley, S 
Programming and Atkinson (1985).  However, I ultimately need to estimate 
these numbers.

	  Thanks,
	  spencer graves