Re: [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
______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html