How to calculate the value of one predictor at which probability of response=0.5 in a multivariate GLMER model? I've developed a glmer model to predict the presence/absence of a certain behavior based on multiple environmental predictors (temperature, windspeed, precipitation). 14 animals were observed continuously throughout the daylight hours over multiple days, and data were summarized over two-hour periods. The response is the proportion of the two-hour period that the individual spent performing the behavior, and the predictors are the average temperature, windspeed, and precipitation over that two-hour period. Because the response variable is a proportion, I've included a weighting factor in the model. I'm most interested in the relationship between temperature and the response behavior. The model is specified as:
Model<-glmer(Response~Temp+Precip+Windspeed+(Temp|Individual), data=data, family=binomial, weights=weightfactor)

I'd like to calculate the temperature at which, with other predictors (precipitation and wind speed) held at their mean values, the probability of the response = 0.5, and I'd also like to calculate 95% confidence intervals for that estimate. I found these threads, which describe how to adapt the dose.p function to work with glmer models. I've used those as a basis for my code:
  dose.p.glmm <-  function(Model, cf = 1:2, p = 0.5) {
  f <- family(Model)
  eta <- f$linkfun(p)
  b <- fixef(Model)[cf]
  x.p <- (eta - b[1L])/b[2L]
  names(x.p) <- paste("p = ", format(p), ":", sep = "")
  pd <- -cbind(1, x.p)/b[2L]
  SE <- sqrt(((pd %*% vcov(Model)[cf, cf]) * pd) %*% c(1, 1))
  res <- structure(x.p, SE = matrix(SE), p = p)
  class(res) <- "glm.dose"
  res}
  dose.p.glmm(Model, cf=1:2, p=0.5)

However, this code was written for glmer models with a single predictor variable. Is there a way to adapt my code to specify that in calculating this, I want to hold precip and windspeed at their mean observed values and then find the value of temp at which the probability of the response=0.5? And if not, is there another way to calculate that value and its 95% confidence intervals?
 A: Your function appears to assume that all other variables are fixed at zero.  (Since they are at zero, neither their coeffs nor their SEs come into the equation).
So, if you would like to calculate the same thing when the other variables are fixed at values other than zero, a quick dirty approach would be to subtract those certain values from the input values for those variables then fit your Model to that recentred data.  Then use your dose.p.glmm function exactly as it is.  Make sure that cf is set to a two element vector that gives the index in fixef(Model) of the intercept and slope that you care about, in that order.  (See the help for MASS::dose.p).
Caveats
A couple of question marks about that approach you should be aware of if you want to use it, however:

*

*It assumes that the values on the "other" variables are fixed values, rather than estimates of means.  It doesn't account for uncertainty in what the means should be.  This may or may not bother you.


*The function was originally written for use in GLMs (in the MASS package).  In GLMMs, we also have random effects, and they're harder to ignore than in linear models, because the nonlinearities mean there are consequences to how we deal with them.  In your adapted function, the random effects are assumed to be zero.  This (I think?) corresponds to the median case (say, if your random effect captures different participants, then we are calculating for the median participant).  But it does not necessarily correspond to the overall probability that you will see -- that would be a mean, not a median, and we would have to integrate over the random effect to get that.  That said, if your p is really going to be 0.5, the symmetry of the logit link around 0.5 should mean that high and low values would cancel each other symmetrically in that integral, so in that specific case, it shouldn't matter either way.
