[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