[R-sig-ME] gee, geese and glmer

[Yang, Qiong]

Hi Ben,

We wonder if you can add an option in lmer() of current lme4 version to call the algorithm used in lme4_0.999999-2?

For our package (analyze rare genetic variant) to be put on CRAN, we need to use current version of lme4. However, at this point, there are still issues that cannot be resolved with newer versions of lme4. It is very difficult for us to keep waiting and testing the new release, and hope all the issues resolved and no new issues coming up. lme4_0.999999-2 has been used by us for a long time with little problem. Your help on this is highly appreciated.
Best,
Qiong

-----Original Message-----
From: Ben Bolker [mailto:bbolker at gmail.com] 
Sent: Saturday, March 08, 2014 5:25 PM
To: r-sig-mixed-models at r-project.org
Cc: Chen, Ming-Huei; Yang, Qiong
Subject: Re: gee, geese and glmer

On 14-03-07 11:25 PM, Ming-Huei Chen wrote:
> Hi Ben,
> 
>  
> 
> In an analysis we found that glmer in new lme4 gave result different 
> from old lme4, gee and geese, where old lme4 seems to be closer to gee 
> and geese.. Please see highlighted sex effect below. Case by sex (2x2) 
> table is also given. Can you please let us know how would you look 
> into the results? Thanks!
> 

   [cc'ing to r-sig-mixed-models: **please** try r-sig-mixed-models first, not personal e-mail to me ...]

  I can't say exactly what's going here; without having a reproducible example <http://tinyurl.com/reproducible-000> it's hard to say precisely.  Thoughts:

 * gee and geese are giving _exactly_ the same parameter estimates, to
8 significant digits, so I would guess they are wrapping identical underlying methods.

 * As far as diagnosing the issue with lme4 1.0-6:
   * does changing the optimization method, i.e.
 glmerControl(optimizer="optimx",optCtrl=list(method="nlminb"))
   [must do library("optimx") first] or
 glmerControl(optimizer="bobyqa")

  change the result?

 * I would be curious whether the soon-to-be-released version 1.1-4 (which can be installed from github or lme4.r-forge.r-project.org/repos) gives either (1) convergence warnings or (2) different/better answers

 * You can try specifying the starting values for lme4 to diagnose misconvergence; for example, start lme4 from the estimates given by old lme4/lme4.0 and see if it gives a similar answer.

 * You can use the 'slice' and 'splom.slice' functions from bbmle to visualize the likelihood surfaces

  good luck,
   Ben Bolker

> Ming-Huei
> 

> ###GEE
> 
>> summary(gee(case~sex+PC1+PC2+PC3+PC4,id=famid,family=binomial,data=da
>> ta))$coef
>                Estimate Naive S.E.     Naive z Robust S.E.   Robust z
> (Intercept) -1.88047373 0.13532162 -13.8963286  0.15960440 -11.782092
> sex         -0.23436854 0.08611269  -2.7216494  0.09050577  -2.589543
> PC1         -0.05478639 0.06195318  -0.8843192  0.06822178  -0.803063
> PC2         -0.09934572 0.06494563  -1.5296753  0.06520811  -1.523518
> PC3         -0.07020391 0.06626875  -1.0593818  0.06962147  -1.008366
> PC4         -0.13413097 0.06746716  -1.9880927  0.06979901  -1.921674
> 

> ###GEESE
> 
>> summary(geese(case~sex+PC1+PC2+PC3+PC4,id=famid,family=binomial,data=
>> data))$mean
> 
>                estimate     san.se        wald           p
> 
> (Intercept) -1.88047373 0.15960440 138.8176912 0.000000000
> sex         -0.23436854 0.09050577   6.7057312 0.009610351
> PC1         -0.05478639 0.06822178   0.6449102 0.421938319
> PC2         -0.09934572 0.06520811   2.3211071 0.127629159
> PC3         -0.07020391 0.06962147   1.0168016 0.313278888
> PC4         -0.13413097 0.06979901   3.6928324 0.054646745
> 
> ### lme4_0.999999-2
> 
>> summary(glmer(case~sex+PC1+PC2+PC3+PC4+(1|famid),family=binomial,data
>> =data))
>             Estimate Std. Error z value Pr(>|z|)
> (Intercept) -3.01599    0.28305 -10.655   <2e-16 ***
> sex         -0.41056    0.16285  -2.521   0.0117 *
> PC1         -0.17116    0.12903  -1.326   0.1847
> PC2         -0.15510    0.13382  -1.159   0.2465
> PC3         -0.19044    0.13580  -1.402   0.1608
> PC4          0.02532    0.13732   0.184   0.8537
> 
> ###lme4_1.0-6
> 
>> summary(glmer(case~sex+PC1+PC2+PC3+PC4+(1|famid),family=binomial,data
>> =data))
> 
>             Estimate Std. Error z value Pr(>|z|)
> 
> (Intercept) -10.2784     0.8631 -11.909   <2e-16 ***
> sex           0.3497     0.1975   1.770   0.0767 .
> PC1          -0.3555     0.1623  -2.190   0.0285 *
> PC2          -0.1087     0.1653  -0.657   0.5109
> PC3          -0.2242     0.1652  -1.357   0.1748
> PC4           0.1103     0.1671   0.660   0.5091
> 
> Case by sex
> 
>         1    2
>   0   2554 3021
>   1    310  290
>