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

Ben Bolker bbolker at gmail.com
Tue Mar 18 15:06:15 CET 2014


  This simplified version should help, I think.
  I'm still puzzled why lme4 1.1-5 is getting stuck with many different
optimizers, and why the slice plot looked sensible.  I would like to be
able to visualize the likelihood surface better.
  Have you gotten bbmle:::slicetrans to work yet?

On 14-03-17 10:49 PM, Ming-Huei Chen wrote:
> Qiong,
> 
> Controls are from extended families but NOT all the cases are singletons.
> Some cases are from extended families. Sorry I didn't make it clear. 
> 
> Best,
> 
> Ming-Huei
> 
> -----Original Message-----
> From: Yang, Qiong [mailto:qyang at bu.edu] 
> Sent: Monday, March 17, 2014 8:18 PM
> To: Ben Bolker; Chen, Ming-Huei
> Cc: 'lme4-authors at lists.r-forge.r-project.org'
> Subject: RE: [R-sig-ME] gee, geese and glmer
> 
> Hi Ben,
> Thanks for be mindful about our data confidentiality- what we communicated
> in this email sequence so far can be copied to the general mailing list.
> 
> Just to clarify further: All the cases are singletons(i.e. families of 1
> member) while controls are from extended families (i.e. families of multiple
> members). 
> 
> The predictor variables do not identify the cases from controls.
> We fitted reduce model where the only predictor is sex which does not
> identify cases from controls however still observed reversed sign on beta
> coefficients between lme4 and other software:
> First, let's look at crude case by sex table
>         1    2
>   0   2554 3021
>   1    310  290
> Ignoring family structure: Odds Ratio=0.79; beta=-0.23
> 
> Since all cases are singletons, I don't think including family-specific
> random effects should change the direction of sex effects. Therefore I tend
> to think the effect direction given by lme4 is not correct. Below are the
> results with only sex as predictor. 
>> library(lme4)
>> m1 <- glmer(case~sex+(1|famid),family=binomial,data=data)
>> m1 at optinfo$warnings
> list()
>> coef(summary(m1))
>                Estimate Std. Error    z value      Pr(>|z|)
> (Intercept) -10.2359223  0.4644516 -22.038725 1.225388e-107
> sex           0.3531757  0.1837722   1.921813  5.462931e-02
>> library(lme4.0)
>> m0 <- lme4.0::glmer(case~sex+(1|famid),family=binomial,data=data)
> Warning message:
> In mer_finalize(ans) : false convergence (8)
>> coef(summary(m0))
>               Estimate Std. Error    z value     Pr(>|z|)
> (Intercept) -3.2807992  0.2868701 -11.436533 2.746424e-30
> sex         -0.2592505  0.1629382  -1.591097 1.115877e-01
> 
>> library(geepack)
>> summary(geese(case~sex,id=famid,family=binomial,data=data))$mean
>               estimate     san.se      wald           p
> (Intercept) -1.8742252 0.16010676 137.03271 0.000000000
> sex         -0.2346185 0.09068416   6.69363 0.009675797
>> library(gee)
>> summary(gee(case~sex,id=famid,family=binomial,data=data))$coef
>               Estimate Naive S.E.    Naive z Robust S.E.   Robust z
> (Intercept) -1.8742252 0.13510959 -13.871889  0.16010676 -11.706097
> sex         -0.2346185 0.08601766  -2.727562  0.09068416  -2.587205
> 
> Thanks a lot, Qiong
> 
> -----Original Message-----
> From: Ben Bolker [mailto:bbolker at gmail.com] 
> Sent: Monday, March 17, 2014 3:51 PM
> To: Chen, Ming-Huei; Yang, Qiong
> Cc: 'lme4-authors at lists.r-forge.r-project.org'
> Subject: Re: [R-sig-ME] gee, geese and glmer
> 
> 
>   [I'm not going to cc: to r-sig-mixed-models, as you may consider parts of
> it to be sensitive, but I would like to ask permission to do so in the
> future, as I think it would be useful to have these conversations in public
> -- please let me know]
> 
> 
> On 14-03-17 03:39 PM, Ming-Huei Chen wrote:
>> Thanks, Ben!
>>
>> Our data was combined by two cohorts, one of them contains only cases 
>> and is an unrelated sample (249 male, 182 female). Do you think this 
>> may lead to what we observed here?
> 
>   Possibly, although I would think the cohorts would have to be identifiable
> somehow from the predictor variables.
> 
>>
>> * I don't know why I got an error...
>> Error: could not find function "slicetrans"
> 
>   slicetrans is a function from the bbmle package, but it turns out to be
> un-exported (this is still rather experimental functionality), so you would
> need to use bbmle:::slicetrans
> 
>>
>> * Do you mean using nAGQ? If so, I used nAGQ=5 an dgot a warning.
>>> x <-
>> summary(lme4::glmer(case~sex+PC1+PC2+PC3+PC4+(1|famid),family=binomial
>> ,data=
>> test.dat,nAGQ=5))
>> Warning message:
>> In checkConv(attr(opt, "derivs"), opt$par, checkCtrl = control$checkConv,
> :
>>   Model failed to converge with max|grad| = 0.00275074 (tol = 0.001) 
>> Fixed effects:
>>              Estimate Std. Error z value Pr(>|z|)
>> (Intercept) -2.675364   0.282841  -9.459   <2e-16 ***
>> sex         -0.053593   0.146444  -0.366   0.7144
>> PC1         -0.205982   0.114760  -1.795   0.0727 .
>> PC2         -0.161422   0.119219  -1.354   0.1757
>> PC3         -0.191577   0.121107  -1.582   0.1137
>> PC4          0.002723   0.122446   0.022   0.9823
> 
>   This warning *might* be spurious (the lme4 maintainers have been having a
> private conversation about the fact that one may get false-positive warnings
> about gradients when the fit is singular):
> 
>>
>> * Do you mean use fixef.prior? fixef.prior: a BLME prior of family 
>> normal, t, or NULL. So which one is weaker prior?
> 
>    t is weaker than Normal (i.e., for the same mean (typically 0) and
> variance it has fatter tails).  However, there are two dimensions of
> "weakness" -- to make a prior weaker you can make the scale (i.e.
> variance) larger, and/or make the shape parameter more diffuse (i.e., make
> the tails fatter, in the sequence Gaussian (= t_\infty -> t_{large n} ->
> t_{smaller n}).  (You may end up discovering that the results depend on the
> strength of the prior ... which will open up interesting cans of worms in
> terms of what the data are actually telling you ...)
>>
>> Best,
>>
>> Ming-Huei
>>
>> -----Original Message-----
>> From: Ben Bolker [mailto:bbolker at gmail.com]
>> Sent: Friday, March 14, 2014 5:33 PM
>> To: Ming-Huei Chen; 'Yang, Qiong'
>> Cc: 'lme4-authors at lists.r-forge.r-project.org'
>> Subject: Re: [R-sig-ME] gee, geese and glmer
>>
>>   OK, this definitely looks well-behaved.  I am leaning still more 
>> heavily toward the idea that new-lme4 is giving an answer that is 
>> technically more correct (the splom plot shows the focal parameters 
>> and the location of the minimum computed 'slice' value as closed and 
>> open circles (I forget which is which); if they don't coincide, that 
>> suggests a misconvergence.  There are a few cases where there is a 
>> barely visible difference (i.e. you can see the open circle peeking 
>> out from around the edge of the closed circle), but they're generally 
>> very close.  The deviance surface looks sensibly smooth/ellipsoidal, and
> none of the parameters appear to be on boundary.
>> (There's a little problem with the slice2D code that limits the ranges 
>> sampled to be either all-positive or all-negative, I think.)
>>
>>    If so, that makes the difference between GEE, Laplace 
>> approximation, and Bayesian GLMM as much of a statistical as a
> computational problem.
>> A few more things I would try to follow up:
>>
>>   * try s2 <- slicetrans(m0, m1, dd); plot(s2) to examine the 
>> cross-section of the deviance function between the new and old parameters
>>   * try using Gauss-Hermite quadrature to see if that makes a difference
>>   * try blme::bglmer to impose a weak prior on the fixed effects (see 
>> if that makes the results closer to MCMCglmm)
>>
>>
>>
>> On 14-03-14 02:15 PM, Ming-Huei Chen wrote:
>>> Here it is.... 
>>>
>>> ###new
>>>> library(lme4)
>>>> m1 <-
>>> glmer(case~sex+PC1+PC2+PC3+PC4+(1|famid),family=binomial,data=test.da
>>> t
>>> )
>>>> m1 at optinfo$warnings
>>> list()
>>>> library(optimx)
>>> Loading required package: numDeriv
>>>> m1.1 <-
>>> glmer(case~sex+PC1+PC2+PC3+PC4+(1|famid),family=binomial,data=test.da
>>> t
>>> ,contr
>>> ol=glmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
>>>> m1.1 at optinfo$warnings
>>> list()
>>>
>>> ###old
>>>> library(lme4.0)
>>>> m0 <-
>>> glmer(case~sex+PC1+PC2+PC3+PC4+(1|famid),family=binomial,data=test.da
>>> t
>>> )
>>> Warning message:
>>> In mer_finalize(ans) : false convergence (8)
>>>> library(optimx)
>>> Loading required package: numDeriv
>>>> m0.0 <-
>>> glmer(case~sex+PC1+PC2+PC3+PC4+(1|famid),family=binomial,data=test.da
>>> t
>>> ,contr
>>> ol=glmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))
>>> Error in do.call(lmerControl, control) :
>>>   could not find function "glmerControl"
>>>
>>> Best,
>>>
>>> Ming-Huei
>>> -----Original Message-----
>>> From: Ben Bolker [mailto:bbolker at gmail.com]
>>> Sent: Friday, March 14, 2014 11:08 AM
>>> To: Ming-Huei Chen; 'Yang, Qiong'
>>> Cc: lme4-authors at lists.r-forge.r-project.org
>>> Subject: Re: [R-sig-ME] gee, geese and glmer
>>>
>>>   The fact that the new deviance is lower suggests that the new 
>>> version of
>>> lme4 is actually finding a _better_ fit to the data; so does the 
>>> existence of a convergence warning from old lme4 (it's interesting 
>>> that using
>>> control=glmerControl(optimizer="optimx",optCtrl=list(method="nlminb")
>>> ) doesn't give you a similar warning ...) . There is still some minor 
>>> possibility that there's a glitch in the deviance calculation -- a 
>>> good cross-check would be to evaluate the deviance of both sets of 
>>> parameters for old lme4 as well, but that's a bit harder to do 
>>> because
>>> old-lme4 doesn't have the same level of flexibility).
>>>
>>>   The parameter comparisons that new-lme4 is giving fairly extreme 
>>> estimates for one of the fixed parameters (old=-3, new=-10; looking 
>>> at the parameter tables below, it looks like this is the intercept) 
>>> and a very large estimate of the random-effects standard deviation. 
>>> (I could have figured this out from the parameter tables too, but I'm 
>>> lazy and graphical presentation makes it much more obvious.)
>>>
>>>   This is looking like a case of (quasi-)complete separation, which 
>>> would explain why MCMCglmm is giving a much less extreme result (it 
>>> imposes weak priors on all of the parameters by default).  It would 
>>> be interesting to try blme::bglmer on this case as well.
>>>
>>> Can you send the (graphical) results of
>>>
>>> dd <- lme4::glmer(...,devFunOnly=TRUE)
>>> library(bbmle)
>>> ss1 <- slice2D(fun = dd, params = m1parms, verbose = FALSE)
>>> splom(ss1)
>>>
>>> ??
>>>
>>> On 14-03-14 09:39 AM, Ming-Huei Chen wrote:
>>>> Thanks, Ben!
>>>>
>>>> The new deviance is lower than the old one. The plot of estimates is 
>>>> also attached.
>>>>> dd(m0parms)
>>>> [1] 3301.035
>>>>> dd(m1parms)
>>>> [1] 2991.881
>>>>
>>>> Best,
>>>>
>>>> Ming-Huei
>>>> -----Original Message-----
>>>> From: Ben Bolker [mailto:bbolker at gmail.com]
>>>> Sent: Thursday, March 13, 2014 11:29 PM
>>>> To: Ming-Huei Chen; 'Yang, Qiong'; 'Martin Maechler'
>>>> Cc: r-sig-mixed-models at r-project.org
>>>> Subject: Re: [R-sig-ME] gee, geese and glmer
>>>>
>>>>   You should be able to install lme4.0 via:
>>>>
>>>> install.packages("lme4.0",
>>>>     repos="http://lme4.r-forge.r-project.org/repos")
>>>>
>>>> It should give identical results to lme4_0.999999-2.
>>>>
>>>>   Also, if you can pull out the coefficients (both the theta and the 
>>>> fixed effect parameters) and save them somewhere, you don't need to 
>>>> run old and new lme4 in the same session.  You just need to be able 
>>>> to pull both sets of parameters into an R session where you have 
>>>> lme4
>>>>> =
>>>> 1.0 loaded, use devFunOnly=TRUE to extract the deviance function, 
>>>> and apply it to both the old and the new sets of parameters.
>>>>
>>>>   Ben Bolker
>>>>
>>>>
>>>> On 14-03-13 08:54 PM, Ming-Huei Chen wrote:
>>>>> Thanks Ben and Martin!
>>>>>
>>>>> I tried lme4 1.1-4 and the result was identical to lme4 1.0-6 
>>>>> without warning (by @optinfo$warnings). As you mentioned about 
>>>>> convergence warning, I redid the analysis with lme4_0.999999-2 and 
>>>>> saw a convergence warning. I also tried MCMCglmm and here is the 
>>>>> result, sex effect is insignificant but the estimate is negative as 
>>>>> in GEE and lme4_0.999999-2. But as you mentioned, it's different
>> methodology..
>>>>>
>>>>> test <-
>>>>> MCMCglmm(case~sex+PC1+PC2+PC3+PC4,random=~famid,family="categorical"
>>>>> ,
>>>>> d
>>>>> ata=da
>>>>> ta)
>>>>>             post.mean  l-95% CI  u-95% CI eff.samp  pMCMC
>>>>> (Intercept) -2.469391 -2.821757 -2.078372    4.290 <0.001 ***
>>>>> sex         -0.055439 -0.330237  0.111188    5.426  0.626
>>>>> PC1         -0.197212 -0.332290 -0.018821    6.833  0.074 .
>>>>> PC2         -0.158433 -0.364921 -0.003554    5.142  0.002 **
>>>>> PC3         -0.130129 -0.320582  0.065973    4.094  0.320
>>>>> PC4          0.002677 -0.306265  0.233621    2.060  0.812
>>>>>
>>>>> As for comparing deviance, because different versions of lme4 need 
>>>>> to be installed in different versions of R and there is actually no 
>>>>> lme4.0, so I don't know how to use two versions of lme4 in the same 
>>>>> R session... Since you planned to get lme4.0 on CRAN, so you should 
>>>>> have the package... in that case, can you please send it to me? Or?
>>>>>
>>>>> Our data contain pedigree information which cannot be shared (even 
>>>>> pedigree structure), so simulate may not help.
>>>>>
>>>>> Best,
>>>>>
>>>>> Ming-Huei
>>>>>
>>>>> -----Original Message-----
>>>>> From: Ben Bolker [mailto:bbolker at gmail.com]
>>>>> Sent: Thursday, March 13, 2014 4:09 PM
>>>>> To: Yang, Qiong; Martin Maechler
>>>>> Cc: r-sig-mixed-models at r-project.org; Chen, Ming-Huei
>>>>> Subject: Re: [R-sig-ME] gee, geese and glmer
>>>>>
>>>>> On 14-03-13 02:01 PM, Yang, Qiong wrote:
>>>>>> Hi Ben and Martin,
>>>>>>
>>>>>> Thanks you for your reply and taking time to explain the situation. 
>>>>>> Maybe I wasn't clear in my previous message or don't have a deeper 
>>>>>> understanding of the CRAN policy: Instead of separate packages for 
>>>>>> "old-CRAN" lme4 (lme4.0), and the current lme4, which is not 
>>>>>> allowed by CRAN maintainer, is it possible to have an option in 
>>>>>> current lme4 to call the algorithm in old lme4? I thought the new
>>>>>> lme4 implements a different algorithm and why not make the old 
>>>>>> algorithm available in current lme4 lmer() as an option. It is 
>>>>>> like we can request pearson or spearman correlations in cor(). 
>>>>>> Sorry if I simplified the problem too much.
>>>>>
>>>>>    It's not so much CRAN policy as the architecture of the software.  
>>>>> A great deal of the internal structure has changed between versions
>>>>> <1.0 and > 1.0, so putting both versions into one package would 
>>>>> basically mean having copies of all of the old and all of the new 
>>>>> code
>>>>> -- and keeping it all safely distinguished and separate would 
>>>>> probably be substantially more work (we have about as much as we 
>>>>> can to do keep up with a single package ...) As we said before, we 
>>>>> were surprised by the rejection of lme4.0 by CRAN -- our previous 
>>>>> planning had assumed that we would be able to put lme4.0 on CRAN.
>>>>> Somewhat to our surprise, there have been relatively few people 
>>>>> contacting us with troubles similar to yours -- for the most part 
>>>>> it seems that is
>>>>> (thankfully) very rare that
>>>>> lme4 gives significantly worse answers than lme4.0.
>>>>>
>>>>>   So far it's not clear that lme4's answers are actually worse than 
>>>>> lme4.0's (I admit that all other things being equal, closer to GEE 
>>>>> is more likely to be correct -- on the other hand, we know that GEE 
>>>>> gives marginal estimates of fixed parameters while GLMM gives 
>>>>> conditional estimates, so we shouldn't expect them to be the same).
>>>>> Can
>>> I confirm:
>>>>>
>>>>> * you've tried this with the most recent development version of 
>>>>> lme4
>>>>> (1.1-4) and you do *not* get any convergence warnings?
>>>>>
>>>>> * have you compared the deviances based on the old (lme4.0 / lme4 <
>>>>> 1.0) and new packages?  Here is the outline of how you would do this:
>>>>>
>>>>> library(lme4)   ## load new package
>>>>> m1 <- glmer(...) ## lme4 fit
>>>>> library(lme4.0)  ## load old package
>>>>> m0 <- lme4.0::glmer(...) ## lme4.0 fit ## extract full parameter 
>>>>> set for lme4.0 and lme4 fits m0parms <-
>>>>> c(getME(m0,"theta"),fixef(m0)) m1parms <-
>>>>> c(lme4::getME(m1,"theta"),fixef(m1))
>>>>> nparms <- length(m0parms)
>>>>> ntheta <- length(getME(m0,"theta"))
>>>>>
>>>>> par(las=1,bty="l")
>>>>> plot(m0parms,m1parms,col=rep(2:1,c(ntheta,nparms-ntheta)))
>>>>> abline(a=0,b=1)
>>>>>
>>>>> ## set up deviance function
>>>>> dd <- lme4::glmer(...,devFunOnly=TRUE) ## compare
>>>>> dd(m0parms)
>>>>> dd(m1parms)
>>>>>
>>>>>   If the new deviance is lower than the old deviance, that's a 
>>>>> strong indication that the new version is actually getting a 
>>>>> _better_ fit than the old one.
>>>>>
>>>>>   Another diagnostic tool is:
>>>>>
>>>>> library(bbmle)
>>>>> ss1 <- slice2D(fun = dd, params = m1parms, verbose = FALSE)
>>>>> splom(ss1)
>>>>>
>>>>>   This will show you cross-sections of the deviance surface ...
>>>>>
>>>>>
>>>>>> To Ben, Sorry that we cannot share data because of confidentiality 
>>>>>> policy. But you are able to debug our problem through PC screen 
>>>>>> sharing where you can gain control and run programs on our Linux 
>>>>>> session without possessing the data. If that sounds like a good 
>>>>>> idea, let us know when you might be available to do so.
>>>>>
>>>>>   Can you generate a simulated example that replicates your problem?
>>>>> simulate() is a useful tool for generating data that are similar, 
>>>>> but not identical to, your original data (you can also anonymize 
>>>>> factor labels etc.).
>>>>>
>>>>>   PC screen sharing seems difficult -- I have easy access to MacOS 
>>>>> and Linux systems, but not Windows.
>>>>>
>>>>>>
>>>>>> Thanks, Qiong -----Original Message----- From: Ben Bolker 
>>>>>> [mailto:bbolker at gmail.com] Sent: Tuesday, March 11, 2014 9:26 AM To:
>>>>>> Martin Maechler; Yang, Qiong Cc: r-sig-mixed-models at r-project.org; 
>>>>>> Chen, Ming-Huei Subject: Re: [R-sig-ME] gee, geese and glmer
>>>>>>
>>>>>> I would also point out that we are indeed very interested, in the 
>>>>>> medium term (the short term is very very busy!), in making sure 
>>>>>> that your issues are resolved.  We would like lme4 to dominate
>>>>>> lme4.0 (i.e., to work better in all circumstances). So far it's 
>>>>>> been a bit difficult since we have been debugging remotely -- 
>>>>>> short of the suggestions I gave below, and without access to a 
>>>>>> reproducible example, it's very hard indeed for me to say much more.
>>>>>>
>>>>>> sincerely Ben Bolker
>>>>>>
>>>>>> On 14-03-11 05:45 AM, Martin Maechler wrote:
>>>>>>>>>>>> Yang, Qiong <qyang at bu.edu> on Mon, 10 Mar 2014 18:19:56
>>>>>>>>>>>> +0000 writes:
>>>>>>>
>>>>>>>> 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?
>>>>>>>
>>>>>>> unfortunately not.
>>>>>>>
>>>>>>> For this reason, we had planned for many months, starting in 
>>>>>>> August
>>>>>>>  2012 and announced on this mailing list at least a year ago that 
>>>>>>> we would provide the 'lme4.0' package  (back-compatible as well 
>>>>>>> as possible in light of computer OS updates incl system 
>>>>>>> libraries, and R updates, ...) in order to provide useRs the 
>>>>>>> possibility of reproducible research and data analysis for their 
>>>>>>> analyses done with "old-CRAN" lme4 versions.
>>>>>>>
>>>>>>> Unfortunately (for us), the CRAN maintainers decided that 
>>>>>>> providing
>>>>>>>  lme4.0 in addition to lme4 was a bad idea, and explicit forbid 
>>>>>>> such actions in the (then new) CRAN policy document.  I'm still 
>>>>>>> not at all happy with that decision.
>>>>>>>
>>>>>>>> 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 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.
>>>>>>>
>>>>>>> Good to hear.  Such cases were exactly the reason why we (lme4
>>>>>>> authors) made such considerable effort to provide  lme4.0 for 
>>>>>>> reproducible research and data analysis.
>>>>>>>
>>>>>>> at the bottom of
>>>>>>> https://github.com/lme4/lme4/blob/master/README.md we mention the 
>>>>>>> state and give installation instruction of lme4.0, but as you 
>>>>>>> say, this does not solve the problem for other package 
>>>>>>> maintainers: If they want a CRAN package, they (currently? I'm 
>>>>>>> optimistic beyond reason :-) cannot have a 'Depends: lme4.0'  (or 
>>>>>>> "Imports:..." or similar).
>>>>>>>
>>>>>>>> Your help on this is highly appreciated.  Best, Qiong
>>>>>>>
>>>>>>> You're welcome; currentl there's not more we can do. Martin
>>>>>>>
>>>>>>> -- Martin Maechler, ETH Zurich, Switzerland
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>> -----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,
>>>>>>>>>> d
>>>>>>>>>> a
>>>>>>>>>> t
>>>>>>>>>> a
>>>>>>>>>> =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=binomia
>>>>>>>>>> l
>>>>>>>>>> ,
>>>>>>>>>> d
>>>>>>>>>> a
>>>>>>>>>> ta=
>>>>>>>
>>>>>>>>>>
>>>>>>>> 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=binomi
>>>>>>>>>> a
>>>>>>>>>> l
>>>>>>>>>> ,
>>>>>>>>>> d
>>>>>>>>>> ata
>>>>>>>
>>>>>>>>>>
>>>>>>>> =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=binomi
>>>>>>>>>> a
>>>>>>>>>> l
>>>>>>>>>> ,
>>>>>>>>>> d
>>>>>>>>>> ata
>>>>>>>
>>>>>>>>>>
>>>>>>>> =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
>>>>>>>>>
>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> R-sig-mixed-models at r-project.org mailing list 
>>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>>>
>>>>>>
>>>>>
>>
> 
>



More information about the R-sig-mixed-models mailing list