Skip to content

Instantly share code, notes, and snippets.

@CerebralMastication
Created September 16, 2010 17:07
Show Gist options
  • Save CerebralMastication/582767 to your computer and use it in GitHub Desktop.
Save CerebralMastication/582767 to your computer and use it in GitHub Desktop.
set.seed(2)
x <- 1:100
y <- 20 + 3 * x
e <- rnorm(100, 0, 60)
y <- 20 + 3 * x + e
plot(x,y)
yx.lm <- lm(y ~ x)
lines(x, predict(yx.lm), col="red")
xy.lm <- lm(x ~ y)
lines(predict(xy.lm), y, col="blue")
# so lm() depends on which variable is x and wich is y
# lm minimizes y distance (the error term is y-yhat)
#normalize means and cbind together
xyNorm <- cbind(x=x-mean(x), y=y-mean(y))
plot(xyNorm)
#covariance
xyCov <- cov(xyNorm)
eigenValues <- eigen(xyCov)$values
eigenVectors <- eigen(xyCov)$vectors
eigenValues
eigenVectors
plot(xyNorm, ylim=c(-200,200), xlim=c(-200,200))
lines(xyNorm[x], eigenVectors[2,1]/eigenVectors[1,1] * xyNorm[x])
lines(xyNorm[x], eigenVectors[2,2]/eigenVectors[1,2] * xyNorm[x])
# the largest eigenValue is the first one
# so that's our principal component.
# but the principal component is in normalized terms (mean=0)
# and we want it back in real terms like our starting data
# so let's denormalize it
plot(x,y)
lines(x, (eigenVectors[2,1]/eigenVectors[1,1] * xyNorm[x]) + mean(y))
# that looks right. line through the middle as expected
# what if we bring back our other two regressions?
lines(x, predict(yx.lm), col="red")
lines(predict(xy.lm), y, col="blue")
@CerebralMastication
Copy link
Author

And Josh Ulrich provided an answer in under 20 minutes. Hive mind FTW!

http://stackoverflow.com/questions/3737165/drop-lines-from-actual-to-modeled-points-in-r/3737183#3737183

-JD

@rubgb
Copy link

rubgb commented Sep 18, 2010

Thanks a lot!

@SwampThingPaul
Copy link

I really enjoyed the blog post and thanks for sharing the code. I have been trying to figure out how to generate the orange (drop) lines on this plot(and below). I tried the segments() function as suggested by Josh Ulrich's posts...but keep going in circles. Any help is appreciated.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment