piecewiseSEM (version 1.2.1)

partial.resid: Calculates partial residuals for two variables

Description

Extracts partial residuals for y ~ x | Z, where Z represents all other variables present in a structured equation upon which x is conditional.

Usage

partial.resid(.formula = y ~ x, modelList, data, model.control = NULL, 
return.data.frame = TRUE, plotit = TRUE, plotreg = TRUE, plotCI = TRUE)

Arguments

.formula

a formula specifying the partial residuals to investigate.

modelList

a single model or list of regressions representing the structural equation model.

data

a data.frame used to construct the model.

model.control

a list of model control arguments to be passed to the partial residual models.

return.data.frame

whether a data.frame of the partial residuals should be returned. Default is TRUE.

plotit

whether the partial plot of y ~ x | Z should be returned. Default is TRUE.

plotreg

whether the partial regression of y ~ x | Z should also be plotted. Default is TRUE.

plotCI

whether the confidence intervals of the partial regression of y ~ x | Z should also be plotted. Default is TRUE.

Value

Returns a data.frame where the first column are the partial residuals of y ~ x | Z, and the second column is the partial residuals of x ~ Z | y.

If plotit = TRUE, then the function also returns a plot of y ~ x | Z.

If plotreg = TRUE, then the plot includes the linear regression of y ~ x | Z.

References

Shipley, Bill. Cause and correlation in biology: a user's guide to path analysis, structural equations and causal inference. Cambridge University Press, 2002.

Examples

# NOT RUN {
# Load model package
library(nlme)

# Load data from Shipley (2013)
data(shipley2013) 

shipley2013.modlist = list(
  
  lme(x2~x1, random = ~x1 | species, data = shipley2013),
  
  lme(x3~x2, random = ~x2 | species, data = shipley2013),
  
  lme(x4~x2, random = ~x2 | species, data = shipley2013),
  
  lme(x5~x3+x4, random = ~x3+x4 | species, data = shipley2013)
  
)

# Get partial residuals of x3 on x5 conditional on x4
resids.df = partial.resid(x5 ~ x3, shipley2013.modlist, shipley2013, 
  list(lmeControl(opt = "optim")))

# Also returns raw residuals values for plotting in other packages
head(resids.df)

# }
# NOT RUN {
  # Create example data
  set.seed(1)
  
  example.data = data.frame(
    y = rnorm(100, 0, 1),
    x1 = rnorm(100, 10, 50),
    random = letters[1:5]
  )
  
  example.data$x2 = example.data$y + runif(100, 0, 5)
  
  example.data$x3 = example.data$x2 + runif(100, 0, 5)
  
  # Run regular linear model using lm()
  lm.model = lm(y ~ x1 + x2 + x3, example.data)
  
  partial.resid(y ~ x3, lm.model, example.data)
  
  # Works with interactions too
  lm.model2 = lm(y ~ x1 * x2 + x3, example.data)
  
  partial.resid(y ~ x1 * x2, lm.model2, example.data, return.data.frame = FALSE)
  
  # Run generalized least squared regression
  library(nlme)
  
  gls.model = gls(y ~ x1 + x2 + x3, example.data)
  
  partial.resid(y ~ x3, gls.model, example.data, return.data.frame = FALSE)
  
  # Run mixed effects model using lme()
  lme.model = lme(y ~ x1 + x2 + x3, random = ~ 1 | random, example.data)
  
  partial.resid(y ~ x3, lme.model, example.data, return.data.frame = FALSE)
  
  # Run mixed effects model using lmer()
  library(lme4)
  
  lmer.model = lmer(y ~ x1 + x2 + x3 + (1 | random), example.data)
  
  partial.resid(y ~ x3, lmer.model, example.data, return.data.frame = FALSE)
  
  # Remove some values of x3 from data.frame
  example.data[c(12, 23, 45), "x3"] = NA
  
  # Run regular linear model using lm()
  lm.model = lm(y ~ x1 + x2 + x3, example.data)
  
  sum(!is.na(partial.resid(y ~ x3, lm.model, example.data)$x.resids)) # Should be 97
  
# }