R: Test for correlation with a covariate? R: I have a linear correlation between two variables x and y, but need to control for a third variable that may be influencing the output. My question is about which test I should/ can use. Can I use a Pearson's correlation and add a covariate? If so, how? Or should I use a GLM with a covariate? How would these approaches differ? I am relatively new to R, so any help on the code would be helpful too. Thank you! 
 A: Yes, you can use correlation with a covariate. This is called partial correlation. It produces a (partial) correlation coefficient that is normalized to the [-1, 1] range just like a regular correlation coefficient, except that the covariate is "controlled for" in the analysis -- a concept which is kind of subtle, but some good explanations of what it really means can be found here.
One way to get a partial correlation in R is using the ppcor package:
# install/load 'ppcor' package for its pcor.test() function
if(!require("ppcor")){
  install.packages("ppcor", repos='http://cran.us.r-project.org')
  library(ppcor)
}

# make up data
x <- rnorm(50)
y <- rnorm(50)
z <- rnorm(50)

# partial correlation between x and y, controlling for z
pcor.test(x, y, z)
#      estimate   p.value  statistic  n gp  Method
# 1 -0.02288511 0.8752972 -0.1569335 50  1 pearson

You also asked whether this differed from using a linear model with a covariate. It doesn't! Check it out:
summary(lm(y ~ x + z))
# Call:
# lm(formula = y ~ x + z)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -2.80457 -0.76631 -0.00539  0.64083  2.79261 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)  0.05163    0.17006   0.304    0.763
# x           -0.02486    0.15842  -0.157    0.876
# z            0.07098    0.15333   0.463    0.646
# 
# Residual standard error: 1.184 on 47 degrees of freedom
# Multiple R-squared:  0.005408,    Adjusted R-squared:  -0.03692 
# F-statistic: 0.1278 on 2 and 47 DF,  p-value: 0.8804

Notice that most of the numbers for the "x" row in the lm() output match those that we got from pcor.test(). The only difference is that the "estimate" for pcor.test() is the partial correlation coefficient, while the "estimate" for lm() is the slope. (The two estimates happen to be numerically similar here, but they are NOT the same.)
