# Am I getting the definition of partial residual plot right?

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
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)