In the context of multiple linear regression, "adjusting" for a covariate simply means including it as an explanatory variable.
There is an equivalent way of understanding multiple linear regression that provides insight into this question. To regress, say, $Z$ on $X$ and $Y$, we may (arbitrarily) select one of the explanatory variables (let it be $X$) and
Separately regress $Z$ on $X$ and $Y$ on $X$, producing residuals $Z_X$ and $Y_X$, respectively, then
Regress $Z_X$ on $Y_X$.
The first step "takes $X$ out of both $Y$ and $Z$." Its interpretation and legitimacy are the same as for any ordinary regression. The generalization to more than two explanatory variables should be clear.
After step 1, you do not have to proceed to step 2: you may study the relationships among the residuals in any ways you please, including PCA. The residuals, by construction, will be orthogonal to a constant (they will have zero means) and also orthogonal to the covariates that have been "left out" ("adjusted for").
You can adjust for interactions and nonlinear terms in the same way: create variables for them and include those in the variables that are taken out.
I will illustrate with R
code.
set.seed(17)
#
# Create independent variables.
#
n.rows <- 30 # Data count
n.vars <- 4 # Number of explanatory variables
n.adj <- 2 # Number of covariates to adjust for
x <- matrix(rnorm(n.rows * n.vars), nrow=n.rows)
#
# Create a dependent variable.
#
intercept <- -1
beta <- 1:n.vars # Coefficients
sigma <- 0.5 # "Error" SD
eps <- rnorm(n.rows, 0, sigma) # "Errors"
y <- intercept + x %*% beta + eps
#
# Multiple regression.
#
#pairs(data.frame(x,y))
summary(lm(y ~ x))
#
# Take out the first `n.adj` variables.
#
x.adj <- x[, 1:n.adj, drop=FALSE] # Variables to take out
x.left <- x[, -(1:n.adj), drop=FALSE] # Variables that will be left
x.res <- apply(x.left, 2, function(y) residuals(lm(y ~ x.adj)))
y.res <- residuals(lm(y ~ x.adj))
#
# Compare the multiple regression of the residuals to the original regression.
#
#pairs(data.frame(x.res, y.res))
summary(lm(y.res ~ x.res))
#
# PCA of the residuals
#
pca <- princomp(cbind(x.res, y.res))
summary(pca)
biplot(pca)
# etc., etc.