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

[Spencer Graves]

Hi, Doug: 

      Thanks.  I'll try the things you suggests.  The observed 
proportions ranged from roughly 0.2 to 0.8 in 100 binomial random 
samples where sigma is at most 0.05.  Jim Lindsey's "glmm" does 
Gauss-Hermite quadrature, but I don't know if it bothers with the 
adaptive step.  With it, I've seen estimates of the variance component 
ranging from 0.4 to 0.7 or so.  Since I simulated normal 0 standard 
deviation of 1, the algorithm was clearly underestimating what was 
simulated.  My next step, I think, is to program adaptive Gauss-Hermite 
quadrature for something closer to my real problem (as you just 
suggested), and see what I get. 

      You mentioned the little vs. big approximations:  My real 
application involves something close to a binomial response driven by 
Poisson defects, where the Poisson defect rate is not constant.  I've 
shown that it can make a difference whether the defect rate is lognormal 
or gamma, so that is another complication and another reason to write my 
own log(likelihood).  I've thought about writing my own function to do 
adaptive Gauss-Hermite quadrature, as you suggested, but decided to 
check more carefully the available tools before I jumped into my own 
software development effort. 

      Thanks again.
      Spencer Graves    

Douglas Bates wrote:

>Spencer Graves <spencer.graves at pdf.com> writes:
>
>  
>
>>	  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
>>    
>>
>
>Check the observed proportions in the data and see if they apparently
>vary enough to be able to expect to estimate a random effect.
>
>It is entirely possible to have the MLE of a variance component be
>zero.
>
>Another thing to do it to check the convergence.  Use
>
>GLMM(y ~ 1, family = binomial, random = ~1|smpl, weigths = n, 
>     control = list(EMv=TRUE, msV=TRUE))
>
>or 
>
>GLMM(y ~ 1, family = binomial, random = ~1|smpl, weigths = n, 
>     control = list(EMv=TRUE, msV=TRUE, opt = 'optim'))
>
>You will see that both optimizers push the precision of the random
>effects to very large values (i.e. the variance going to zero) in the
>second of the penalized least squares steps.
>
>I think that this is a legitimate optimum for the approximate
>problem.  It may be an indication that the approximate problem is not
>the best one to use.  As George Box would tell us,
>
> You have a big approximation and a small approximation.  The big
> approximation is that your approximation to the problem you want to
> solve.  The small approximation is involved in getting the solution
> to the approximate problem.
>
>For this case, even if I turn off the PQL iterations and go directly
>to the Laplacian approximation I still get a near-zero estimate of the
>variance component.  You can see the gory details with
>
>GLMM(y ~ 1, family = binomial, random = ~1|smpl, weigths = n,
>     control = list(EMv=TRUE, msV=TRUE, glmmMaxIter = 1), method = 'Laplace')
>
>I am beginning to suspect that for these data the MLE of the variance
>component is zero.
>
>  
>
>>	  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.
>>    
>>
>
>I think the most reliable method of fitting this particular form of a
>GLMM is using adaptive Gauss-Hermite quadrature to evaluate the
>marginal likelihood and optimizing that directly.  In this model the
>marginal likelihood is a function of two parameters.  If you have
>access to SAS I would try to fit these data with PROC NLMIXED and see
>what that does.  You may also be able to use Goran Brostrom's package
>for R on this model.  As a third option you could set up evaluation of
>the marginal likelihood using either the Laplacian approximation to
>the integral or your own version of adaptive Gauss-Hermite and look at
>the contours of the marginal log-likelihood.
>
>______________________________________________
>R-help at stat.math.ethz.ch 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
>  
>