I am trying to make a partial residual plot of one quadratic predictor in GLMM (
glmer in R). I tried to do it manually, from how I understood the definition of partial residual plot, but I'm still getting different result than the one of the
effects package. Where do I get it wrong?
This is my GLMM model:
> m3b <- glmer(cbind(Juv, Ad) ~ cov1 + I(cov1^2) + Ad_scaled + (1|Species:Year) + (1 + cov1 + I(cov1^2)|Species) + (1|Species:Site), family = binomial, data = data, control = glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb'))) > summary(m3b) Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod'] Family: binomial ( logit ) Formula: cbind(Juv, Ad) ~ cov1 + I(cov1^2) + Ad_scaled + (1 | Species:Year) + (1 + cov1 + I(cov1^2) | Species) + (1 | Species:Site) Data: data Control: glmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb")) AIC BIC logLik deviance df.resid 2970.6 3026.3 -1473.3 2946.6 754 Scaled residuals: Min 1Q Median 3Q Max -4.2727 -0.8569 -0.2288 0.5996 3.2889 Random effects: Groups Name Variance Std.Dev. Corr Species:Site (Intercept) 0.475717 0.68972 Species:Year (Intercept) 0.023712 0.15399 Species (Intercept) 0.980009 0.98995 cov1 0.002376 0.04874 1.00 I(cov1^2) 0.003590 0.05992 -1.00 -1.00 Number of obs: 766, groups: Species:Site, 160; Species:Year, 40; Species, 4 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.08008 0.50172 0.160 0.8732 cov1 -0.13683 0.04656 -2.939 0.0033 ** I(cov1^2) 0.11609 0.05035 2.306 0.0211 * Ad_scaled -0.29382 0.02932 -10.022 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects: (Intr) cov1 I(1^2) cov1 0.533 I(cov1^2) -0.645 -0.561 Ad_scaled -0.007 0.031 -0.108 convergence code: 0 Parameters or bounds appear to have different scalings. This can cause poor performance in optimization. It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
So now, as I understand the definition of partial residual plot, I make my own code according to the definition:
plot(data$cov1, residuals(m3b, type = "response") + data$cov1*-0.13683 + data$cov1^2*0.11609) curve(x*-0.13683 + x^2*0.11609, xlim = c(min(data$cov1), max(data$cov1)), add = TRUE)
But when I use prepared functions in package
effects, the plot is different! The quadratic curve looks similar in both (apart from scaling, which is OK), but there are apparent and significant differences in the individual points. How is it possible? I guess the
effects package has it right. Where did I strayed from the correct definition of partial plot?
require(effects) est<-Effect("cov1", partial.residuals=TRUE, m3b) plot(est)