On 10-09-09 06:03 PM, Russell Millar wrote:
> Dear All,
>
> I'm fitting to some historical grouped binomial data. There are 16
> binomial obs arranged in 8 blocks of 2 (trmt and control). I know
> what the model deviance calculated via the Laplace approx and
> quadrature should be since I've fitted the model using SAS PROC
> GLIMMIX, ADMB and via simulated likelihood (Millar 2004, ANZJS).
> However, lmer is only giving the correct deviance under the
> Laplace approximation.  Also, the deviance calculated by lmer
> using quadrature depends on the choice of quadrature points???
>
> By the way, I'm not directly interested in the deviance. The
> calculated value of 35.8 on 13 dof would suggest lack of fit,
> which is nonsense. I've done a parametric bootstrap and this
> observed deviance is typical. I'm really interested in the accuracy
> of the calculated log-likelihood.
>
> Thanks for any help,
>
> Russell Millar Dept of Stats U. Akld
>
>
> #Binomial data in 8 blocks of 2 clinic=rep(1:8,rep(2,8))
> trmt=as.factor(rep(c("Drug","Control"),8))
> y=c(11,10,16,22,14,7,2,1,6,0,1,0,1,1,4,6) #Successes
> n=c(36,37,20,32,19,19,16,17,17,12,11,10,5,9,6,7) #Number of trials
>
> #log-likelihood of saturated model sum(log(dbinom(y,n,y/n))) #[1]
> -19.19493
>
> #Using 3 forms of alternative software, it has been calculated
> that #the true log-likelihood (with consts) of this model is
> -37.0313 #or -37.1006 if calculated using the Laplace
> approximation. #These correspond to deviances of
> #2*(-19.19493+37.0313)=35.67274 (exact)
> #2*(-19.19493+37.1006)=35.81134 (from Laplace approx)
>

  Confirmed.  Just for bloody-mindedness, I calculated deviances for
every value of
nAGQ from 2 to 99.  (I tried glmmML too -- it gives 35.81135 for all
choices of number of points, and is actually much faster than lmer for
this particular problem ...)

sfun <- function(npts) {
  deviance(
lmer(y/n~trmt+(1|clinic),family="binomial",weight=n,nAGQ=npts) )
}

ss <- sapply(2:99,sfun)

Things get wacky beyond nAGQ=80 ...

 dput(ss)
c(38.9024410477493, 38.9024410477493, 39.0462048218837, 39.0462048218837,
39.0513007237354, 39.0513007237354, 39.0513946823561, 39.0513946823561,
39.0513897918250, 39.0513897918250, 39.0513892605341, 39.0513892605341,
39.0513892516897, 39.0513892516897, 39.0513892540625, 39.0513892540625,
39.0513892542464, 39.0513892542464, 39.0513892542331, 39.0513892542331,
39.0513892542298, 39.0513892542298, 39.0513892542296, 39.0513892542296,
39.0513892542296, 39.0513892542296, 39.0513892542298, 39.0513892542298,
39.0513892542297, 39.0513892542297, 39.0513892542298, 39.0513892542298,
39.0513892542297, 39.0513892542297, 39.0513892542297, 39.0513892542297,
39.0513892542297, 39.0513892542297, 39.0513892542298, 39.0513892542298,
39.0513892542297, 39.0513892542297, 39.0513892542296, 39.0513892542296,
39.0513892542298, 39.0513892542298, 39.0513892542297, 39.0513892542297,
39.0513892542298, 39.0513892542298, 39.0513892542297, 39.0513892542297,
39.0513892542297, 39.0513892542297, 39.0513892542298, 39.0513892542298,
39.0513892542298, 39.0513892542298, 39.0513892542297, 39.0513892542297,
39.0513892542298, 39.0513892542298, 39.0513892542297, 39.0513892542297,
39.0513892542297, 39.0513892542297, 39.0513892542297, 39.0513892542297,
39.0513892542297, 39.0513892542297, 39.0513892542296, 39.0513892542296,
39.0513892542297, 39.0513892542297, 39.0513892542296, 39.0513892542296,
39.0513892542297, 39.0513892542297, 45.0003711126431, 45.0003711126431,
67.5976069661895, 67.5976069661895, 44.860265377668, 44.860265377668,
26.6996251097477, 26.6996251097477, 39.0513892542298, 39.0513892542298,
26.9194621454332, 26.9194621454332, 31.4246100115320, 31.4246100115320,
26.7573461186564, 26.7573461186564, 26.7639453535308, 26.7639453535308,
31.5959347317676, 31.5959347317676)

  A little bit more poking around shows that it's not just the
deviance calculation --
the fixed-effect parameter estimates (at least) vary in a parallel
way, although only in the fourth decimal
place or so.



	[[alternative HTML version deleted]]


