[R-sig-ME] Fwd: Re: Plotting partial residuals from a glmmADMB model

[Ben Bolker]

Here are some examples of partial residual plots ... the only fanciness
available in the base-R machinery [e.g. for partial residuals for lm()
objects] that can't be implemented with the simple machinery here is
collapsing entire terms (due to e.g. factors with more than two levels,
orthogonal polynomials, splines ...) into a single effect

library(glmmADMB)
## ran without FoodTreatment, had trouble with convergence ...
owlfit <- glmmadmb(SiblingNegotiation~ArrivalTime+SexParent+
                       (1|Nest),family="nbinom1",data=Owls)

X <- model.matrix(~ArrivalTime+SexParent,data=Owls)

beta <- fixef(owlfit)

## <https://en.wikipedia.org/wiki/Partial_residual_plot> gives the
## definition of a partial residual as (residuals + \hat beta_i X_i)

## multiply each column of X by the corresponding beta
beta_X <- sweep(X,MARGIN=2,STATS=beta,FUN="*")

## add residuals (columnwise addition works automatically,
## although we could use sweep with MARGIN=1)
p_resid <- sweep(beta_X,MARGIN=1,STATS=residuals(owlfit),FUN="+")

par(mfrow=c(1,2))
for (i in 2:3) {
    plot(X[,i],p_resid[,i],xlab=colnames(X)[i],ylab="partial residuals")
}

(Would be better to plot partial residuals for the binary predictor
as a boxplot ...)

If you were doing this in glmmTMB you can take a couple of shortcuts:

library(glmmTMB)
owlfit2 <- glmmTMB(SiblingNegotiation~ArrivalTime+SexParent+FoodTreatment+
                       (1|Nest),family="nbinom1",data=Owls)

X <- getME(owlfit2,"X")
beta <- fixef(owlfit2)$cond
beta_X <- sweep(X,MARGIN=2,STATS=beta,FUN="*")
p_resid <- sweep(beta_X,MARGIN=1,STATS=residuals(owlfit2),FUN="+")
par(mfrow=c(1,3))
for (i in 2:4) {
    plot(X[,i],p_resid[,i],xlab=colnames(X)[i],ylab="partial residuals")
}



> -----Original Message----- From: Sivakoff, Frances Sent: Tuesday,
> December 19, 2017 3:42 PM To: 'Ben Bolker' <bbolker at gmail.com> Cc:
> r-sig-mixed-models at r-project.org Subject: RE: [R-sig-ME] Plotting
> partial residuals from a glmmADMB model
> 
> Dear Dr. Bolker, Thank you very much for your response. After trying
> to implement your suggestion, I'm unfortunately stuck. It appears
> that the getME function does not work with glmmADMB. I get the
> following error message:
> 
> Error in UseMethod("getME") : no applicable method for 'getME'
> applied to an object of class "glmmadmb"
> 
> A potential work around to this may be to use the "predict" function,
> which can generate the components, but I'm not sure if this is
> equivalent. Also, I'm having trouble following the steps that you
> outlined in your email to generate the partial residuals. Would you
> be willing to work through how to generate the partial residuals for
> the fixed effect "FoodTreatment" in the model below that uses the Owl
> data set?
> 
> ##Using the Owl Data om <-
> glmmadmb(SiblingNegotiation~FoodTreatment+ArrivalTime+SexParent+
> (1|Nest),family="nbinom1",data=Owls)
> 
> Thank you, Frances
> 
> -----Original Message----- From: Ben Bolker
> [mailto:bbolker at gmail.com] Sent: Thursday, December 14, 2017 9:54 PM 
> To: Sivakoff, Frances <sivakoff.3 at osu.edu> Cc:
> r-sig-mixed-models at r-project.org Subject: Re: [R-sig-ME] Plotting
> partial residuals from a glmmADMB model
> 
> You'll have to find a way to make "partial predictions". I don't
> think there's anything built in for this.  *Very* briefly,
> considering only fixed effects, if you retrieve the X matrix 
> (getME(fitted_model,"X")) and the fixed-effect parameters
> (fixef(fitted_model)), you can drop any columns/parameters you want
> and do exp(X %*% beta) with the remaining columns/parameters to get a
> prediction that includes some but not all of the predictors.
> Subtracting the observed value should get you the partial residuals
> ...
> 
> 
> On Tue, Dec 12, 2017 at 2:19 PM, Sivakoff, Frances
> <sivakoff.3 at osu.edu> wrote:
>> I would like to use ggplot2 to plot the partial residuals of an
>> indicator (0 or 1) independent variable in a generalized linear
>> mixed model fit with a "nbinom1" family using glmmadmb. My model
>> has a response variable that is a count, 3 explanatory variables
>> that are continuous, and 7 indicator variables that are 0 when a
>> particular heavy metal is not detected and 1 when it is detected
>> above a threshold value. I'd like to plot the partial residuals of
>> the various independent variables. I think that the model below
>> using the Owl data would be a good example data set for how to do
>> this. The model below has a response variable that is a count, an
>> explanatory variables that is continuous (arrivalTime), two
>> categorical variables (FoodTreatment and SexParent), a random
>> effect of Nest, and uses a "nbinom1" family.
>> 
>> ##Using the Owl Data om <-
>> glmmadmb(SiblingNegotiation~FoodTreatment+ArrivalTime+SexParent+ 
>> (1|Nest),family="nbinom1",data=Owls)
>> 
>> Could you please suggest a method for plotting the partial
>> residuals of the explanatory variables.
>> 
>> Thank you, Frances
>> 
>> 
>> [[alternative HTML version deleted]]
>> 
>> _______________________________________________ 
>> R-sig-mixed-models at r-project.org mailing list 
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models