# Exact synthetic illustration; base R only. Orthogonal residuals.
x <- rep(c(-1,1),each=4)
w <- rep(c(-1,-1,1,1),2)
e <- rep(c(-1,1),4)
z <- x+w
y <- -x+2*z+e
d <- data.frame(x,z,y)
a <- lm(y~x,d); b <- lm(y~x+z,d)
cat(sprintf('Same rows: unadjusted n=%d; adjusted n=%d\n',nobs(a),nobs(b)))
cat(sprintf('Unadjusted x=%.3f; adjusted x=%.3f; adjusted z=%.3f\n',coef(a)['x'],coef(b)['x'],coef(b)['z']))
vif <- 1/(1-summary(lm(x~z,d))$r.squared)
cat(sprintf('Correlation x,z=%.3f; x VIF=%.3f\n',cor(x,z),vif))
cat(sprintf('Unadjusted x CI=[%.3f,%.3f]; adjusted x CI=[%.3f,%.3f]\n',confint(a)['x',1],confint(a)['x',2],confint(b)['x',1],confint(b)['x',2]))
stopifnot(abs(coef(a)['x']-1)<1e-8,abs(coef(b)['x']+1)<1e-8)