Estimating coefficients of an interaction in a GLMM I fitted a Generalized Linear Mixed Model with a Gamma distribution of errors on lme4. I included an interaction between two continuous variables (x:z). How can I estimate the coefficients (beta, Standard Error) of the first variable (x) for different values of the second (z)? I would like to report changes in the slope of x in a paper.
 A: This is really a general question about linear models. Suppose that your linear predictor is
$$
\eta = \beta_0 + \beta_1 x + \beta_2 y + \beta_3 xz
$$
that is, $\beta_3$ is the interaction parameter.  I think what you mean by "the coefficient of $x$ for different values of $z$" is $\beta_1 + \beta_3 z$, i.e. rewriting the expression as 
$$
\beta_0 + (\beta_1 + \beta_3 z) x + \beta_2 z .
$$
In other words, $\beta_1 + \beta_3 z$ gives the expected change in the response (on the linear predictor scale) caused by a 1-unit increase in $x$ at a particular value of $z$.
If you know the variances of $\beta_1$ and $\beta_3$ and their covariance, you can also calculate the variance (and hence standard deviation) of this value. Suppose you have a contrast vector: in this case it would be $c=(0~1~0~z)$, since you're interested in the expression 
$$
\beta_1 + \beta_3 z = 0\cdot \beta_0 + 1 \cdot \beta_1 + 0 \cdot \beta_2 + z \cdot \beta_3
$$
then the variance of the expression is $c \mathbf V c^\top$.
Given that you can specify a focal value of z you're interested in, in R code, you would use
z <- ?? ## a focal value
coefs <- fixef(model)
cc <- coefs["x"] + coefs["x:z"]*z
varcov <- vcov(model)
vv <- varcov[c("x","x:z"),c("x","x:z")]
res_var <- c(1,y) %*% vv %*% c(1,z)  ## variance of the result
res_sd <- sqrt(res_var)

