3
$\begingroup$

I have a reciprocal transplant experiment with larvae from four locations and their substrates with three replicates (blocks). A significant interaction between location and substrate would suggest local adaptation.

I would like to construct a mixed model to test two things: 1) if survival is greater in 'locals' versus foreigners, and 2) if survival is greater for genotypes in their "home" locations versus "away" (substrates). I would use Tukey posthoc tests to determine the above two criteria for local adaptation.

Sample dataset:

structure(list(Location = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 3L, 3L, 3L, 
3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 3L), .Label = c("A", "B", 
"C", "D"), class = "factor"), Substrate = structure(c(2L, 2L, 
2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 4L, 4L, 4L, 
3L, 3L, 3L, 4L, 4L, 4L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 
1L, 3L, 3L, 3L, 4L, 4L, 3L, 3L, 3L, 4L, 4L, 4L, 2L, 4L), .Label = c("A", 
"B", "C", "D"), class = "factor"), Replicate = c(1, 2, 3, 1, 
2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 
2, 3, 1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 3, 2, 3, 1, 2, 3, 1, 2, 3, 
3, 1), Survival = structure(c(5, 10, 14, 18, 27, 18, 1, 22, 9, 
30, 18, 30, 25, 27, 18, 19, 21, 19, 19, 23, 23, 4, 17, 19, 9, 
15, 9, 19, 27, 8, 9, 7, 13, 29, 26, 8, 22, 18, 15, 18, 19, 29, 
15, 30, 30, 30, 0, 0, 25, 20, 16, 12, 3, 12, 29, 8, 21, 0, 12, 
0, 5, 3, 12, 11, 9, 11, 11, 7, 7, 26, 13, 11, 21, 15, 21, 11, 
3, 22, 21, 23, 17, 1, 4, 22, 8, 12, 15, 12, 11, 1, 15, 0, 0, 
0, 30, 30), .Dim = c(48L, 2L))), .Names = c("Location", "Substrate", 
"Replicate", "Survival"), row.names = c(NA, 48L), class = "data.frame")

I have constructed a mixed model as follows:

require(lme4)
Survival.model.1 <- glmer(Survival ~ Location + Substrate + Location:Substrate  + (1|Replicate), data = model_data, family=binomial)

When I run the model I get the following:

Warning messages: 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :   unable to evaluate scaled gradient 
                  2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

The summary of the model results shows that one of the effects has a very large standard error:

Fixed effects:
                      Estimate Std. Error z value Pr(>|z|)
(Intercept)             1.9591     0.4257   4.602 4.19e-06 ***
LocationB              -1.0619     0.3926  -2.705  0.00684 **
LocationC              -1.5271     0.3844  -3.973 7.09e-05 ***
LocationD              -0.7680     0.4017  -1.912  0.05586 .
SubstrateB             -2.5875     0.3882  -6.665 2.65e-11 ***
SubstrateC             -0.9485     0.3957  -2.397  0.01653 *
SubstrateD             -2.1938     0.3832  -5.724 1.04e-08 ***
LocationB:SubstrateB    0.9046     0.5088   1.778  0.07544 .
LocationC:SubstrateB    1.5780     0.4995   3.159  0.00158 **
LocationD:SubstrateB   -0.2139     0.5399  -0.396  0.69190
LocationB:SubstrateC    1.3718     0.5280   2.598  0.00938 **
LocationC:SubstrateC    0.6604     0.5026   1.314  0.18885
LocationD:SubstrateC    0.6546     0.5240   1.249  0.21157
LocationB:SubstrateD    1.9798     0.5039   3.929 8.52e-05 ***
LocationC:SubstrateD    1.1843     0.4958   2.389  0.01691 *
LocationD:SubstrateD   21.5690  2915.9603   0.007  0.99410
--- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation matrix not shown by default, as p = 16 > 12. Use print(x, correlation=TRUE)  or      vcov(x)     if you need it
convergence code: 0 unable to evaluate scaled gradient Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

I believe this is due to the fact that for the three observations of the LocationD x SubstrateD there is 100% survival, while this does not occur anywhere else in our study.

I have found a couple of similar examples here and here an excellent explanation with examples by Ben Bolker here of complete separation or maybe it is called Hauck-Donner effect.

If I follow the example Ben Bolker suggests using:

require(blme)
Survival.model.2 <- bglmer(Survival ~ Location + Substrate + Location:Substrate  + (1|Replicate), data = model_data, family=binomial, fixef.prior = normal(cov = diag(X,16)))

Due to my relatively new experiences with mixed-models, I am unsure how to estimate the prior variance of each parameter (X) in the above code. Or alternatively, are there reasonable default boundary values that can be provided for these parameters?

Any suggestions on how to solve my problem with complete separation is very much appreciated.

$\endgroup$
1
$\begingroup$

With help from others, I have figured out the estimate for prior variance.

Survival.model.2 <- bglmer(Survival ~ -1 + Location:Substrate + (1|Replicate), model_data, family=binomial, fixef.prior = normal(cov = diag(7,16)), control=glmerControl(optimizer="bobyqa"))

In this case though I am still left with the question of whether the random effect is truly important to the model and whether the more simple approach would be to use a Firth logistic regression model?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.