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.
