gee, geese and glmer

Thanks Ben!

S2 plot is attached. 

I got an error using blme...

blme::bglmer(case~sex+PC1+PC2+PC3+PC4+(1|famid),family=binomial,data=test.da
t,fixef.prior=t)
Error in get("checkConv", lme4Env)(attr(opt, "derivs"), opt$par, ctrl =
control$checkConv,  :
  unused argument (ctrl = control$checkConv)

Best,

Ming-Huei

-----Original Message-----
From: Ben Bolker [mailto:bbolker at gmail.com] 
Sent: Monday, March 17, 2014 3:51 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


  [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