Correlation demonstration!

Let’s approach the correlation r again, but from a slightly different angle. I want to demonstrate a helpful way to visualize the calculation.

As before, divide the scatterplot into quadrants (whose axes intersect at the mean of each variable, in this case 0). Now, for each point (i.e., for each pair of variables) find how far it is away from the x-axis and the y-axis.

data<-data.frame(scale(data))
knitr::kable(round(data,3),align='cc')

p<-ggplot(data,aes(X,Y))+geom_point()+geom_hline(yintercept=0)+geom_vline(xintercept=0)+scale_x_continuous(limits=c(-2.5,2.5))+scale_y_continuous(limits = c(-2.5,2.5))+xlab("Zx")+ylab("Zy")
p+geom_segment(aes(x=data[9,1],xend=data[9,1],y=0,yend=data[9,2]),lty=3)+geom_segment(aes(x=0,xend=data[9,1],y=data[9,2],yend=data[9,2]),lty=3)+xlab("Zx")+ylab("Zy")

p+geom_segment(aes(x=data[9,1],xend=data[9,1],y=0,yend=data[9,2]),lty=3)+geom_segment(aes(x=0,xend=data[9,1],y=data[9,2],yend=data[9,2]),lty=3)+annotate("text", x = data[9,1], y=data[9,2], label = paste("(",round(data[9,1],2),",",round(data[9,2],2),")",sep=""),hjust=0,vjust=0,size=3)+xlab("Zx")+ylab("Zy")

Now let’s go ahead and do this for every point:

p+geom_text(aes(label=paste("(",round(data$X,2),",",round(data$Y,2),")",sep=""),hjust=0, vjust=0),size=3)+xlab("Zx")+ylab("Zy")

Since the mean of Zx is ˉx=0 and the mean of Zy is ˉy=0, for each point (xi,yi) we just found the deviations of each value from the mean (xi−ˉx,yi−ˉy)!

Now then, for each pair of values, multiply the two numbers together. Notice that this is the deviation of xi from its mean times the deviation of yi from its mean for all of the i pairs in our sample.

p+geom_text(aes(label=paste("(",round(data$X*data$Y,2),")",sep=""),hjust=0, vjust=0),size=3)+xlab("Zx")+ylab("Zy")

Let’s add them up in each quadrant and put the total in the middle:

p1<-ggplot(data,aes(X,Y))+geom_abline(slope = 1,intercept=0,color='gray50',lty=2)+geom_point()+geom_hline(yintercept=0)+geom_vline(xintercept=0)+scale_x_continuous(limits=c(-3,3))+scale_y_continuous(limits = c(-3,3))+xlab("Zx")+ylab("Zy")+annotate("text",2,2,label=round(sum(data[data$X>0 & data$Y>0,]$X*data[data$X>0 & data$Y>0,]$Y),2))+
  annotate("text",-2,-2,label=round(sum(data[data$X<0 & data$Y<0,]$X*data[data$X<0 & data$Y<0,]$Y),2))+
  annotate("text",-2,2,label=round(sum(data[data$X<0 & data$Y>0,]$X*data[data$X<0 & data$Y>0,]$Y),2))+
  annotate("text",2,-2,label=round(sum(data[data$X>0 & data$Y<0,]$X*data[data$X>0 & data$Y<0,]$Y),2))+xlab("Zx")+ylab("Zy")+geom_rug(color="gray50")
p1

Now let’s use these totals to calculate the correlation coefficient: add up the total in each quadrant and divide by n−1.

∑ZxZyn−1=4+4+0+08=1

Perfect correlation (r=1)! Also, notice that the slope of the (dashed) line that the points lie on is also exactly 1…

What if the relationship were flipped?

data1<-data
data1$X<-data$X*-1
ggplot(data1,aes(X,Y))+geom_point()+geom_hline(yintercept=0)+geom_vline(xintercept=0)+scale_x_continuous(limits=c(-3,3))+scale_y_continuous(limits = c(-3,3))+xlab("Zx")+ylab("Zy")+annotate("text",2,2,label=round(sum(data1[data1$X>0 & data1$Y>0,]$X*data1[data1$X>0 & data1$Y>0,]$Y),2))+
  annotate("text",-2,-2,label=round(sum(data1[data1$X<0 & data1$Y<0,]$X*data1[data1$X<0 & data1$Y<0,]$Y),2))+
  annotate("text",-2,2,label=round(sum(data1[data1$X<0 & data1$Y>0,]$X*data1[data1$X<0 & data1$Y>0,]$Y),2))+
  annotate("text",2,-2,label=round(sum(data1[data1$X>0 & data1$Y<0,]$X*data1[data1$X>0 & data1$Y<0,]$Y),2))+xlab("Zx")+ylab("Zy")+geom_rug(color="gray50")

The sum of the deviations are the same, but the sign has changed. Thus, the correlation is the same strength, but different direction (r=-1) We have a perfect negative correlation: r=∑ZxZyn−1=−88=−1. Notice that the mean and standard deviation of x and y are still the same (we just rotated it 90o)

Now let’s see what would happen if data points were no longer perfectly linear. Let’s take the second data point from the bottom and scoot it up, making the y-value equal 0.

data1<-data
data1[2,2]<-0
ggplot(data1,aes(X,Y))+geom_point()+geom_hline(yintercept=0)+geom_vline(xintercept=0)+scale_x_continuous(limits=c(-3,3))+scale_y_continuous(limits = c(-3,3))+xlab("Zx")+ylab("y")+annotate("text",2,2,label=round(sum(data1[data1$X>0 & data1$Y>0,]$X*data1[data1$X>0 & data1$Y>0,]$Y),2))+
  annotate("text",-2,-2,label=round(sum(data1[data1$X<0 & data1$Y<0,]$X*data1[data1$X<0 & data1$Y<0,]$Y),2))+
  annotate("text",-2,2,label=round(sum(data1[data1$X<0 & data1$Y>0,]$X*data1[data1$X<0 & data1$Y>0,]$Y),2))+
  annotate("text",2,-2,label=round(sum(data1[data1$X>0 & data1$Y<0,]$X*data1[data1$X>0 & data1$Y<0,]$Y),2))+geom_rug(color="gray50")

ggplot(data1,aes(X,Y))+geom_abline(slope = .92,intercept=0,color='gray50',lty=3)+
  geom_abline(slope = 1,intercept=0,color='gray50',lty=2)+geom_point()+geom_hline(yintercept=0)+geom_vline(xintercept=0)+scale_x_continuous(limits=c(-3,3))+scale_y_continuous(limits = c(-3,3))+xlab("Zx")+ylab("Zy")+annotate("text",2,2,label=round(sum(data1[data1$X>0 & data1$Y>0,]$X*data1[data1$X>0 & data1$Y>0,]$Y),2))+
  annotate("text",-2,-2,label=round(sum(data1[data1$X<0 & data1$Y<0,]$X*data1[data1$X<0 & data1$Y<0,]$Y),2))+
  annotate("text",-2,2,label=round(sum(data1[data1$X<0 & data1$Y>0,]$X*data1[data1$X<0 & data1$Y>0,]$Y),2))+
  annotate("text",2,-2,label=round(sum(data1[data1$X>0 & data1$Y<0,]$X*data1[data1$X>0 & data1$Y<0,]$Y),2))+xlab("Zx")+ylab("Zy")+geom_rug(color="gray50")

Uh oh, but since we changed a y-value, the mean of the ys no longer equals 1, and the standard deviation no longer equals 0. They are now My=.12 and SDy=.91; however, X still remains standardized because nothing changed (see distribution along the bottom).

If we ignore this fact and try to calculate the “correlation” we get: rfake=∑ZxYn−1=6.88=.85. However, if we go to the data and calculate the actual correlation, we get this:

knitr::kable(round(data1,3),align='cc')

print(paste("Correlation of X and Y: ",with(data1,cor(X,Y))))

## [1] "Correlation of X and Y:  0.931128347758783"

So the actual r=.93. Moving one data point up reduced our actual correlation from 1 to .93. But what is the relationship between the “fake” correlation we got by not standardizing Y (rfake=.85), and the actual correlation (.93)?

rfake=r×x.85=.93xx=.93.85=.91=sdY

Our “fake” correlation is just the covariance of X and Y! Remember, r=cov(X,Y)sdXsdY

Let’s standardize Y again and proceed:

data1$Y<-scale(data1$Y)
ggplot(data1,aes(X,Y))+geom_abline(slope = .9311,intercept=0,color='gray50',lty=3)+
  geom_abline(slope = 1,intercept=0,color='gray50',lty=2)+geom_point()+geom_hline(yintercept=0)+geom_vline(xintercept=0)+scale_x_continuous(limits=c(-3,3))+scale_y_continuous(limits = c(-3,3))+xlab("Zx")+ylab("Zy")+annotate("text",2,2,label=round(sum(data1[data1$X>0 & data1$Y>0,]$X*data1[data1$X>0 & data1$Y>0,]$Y),2))+
  annotate("text",-2,-2,label=round(sum(data1[data1$X<0 & data1$Y<0,]$X*data1[data1$X<0 & data1$Y<0,]$Y),2))+
  annotate("text",-2,2,label=round(sum(data1[data1$X<0 & data1$Y>0,]$X*data1[data1$X<0 & data1$Y>0,]$Y),2))+
  annotate("text",2,-2,label=round(sum(data1[data1$X>0 & data1$Y<0,]$X*data1[data1$X>0 & data1$Y<0,]$Y),2))+xlab("Zx")+ylab("Zy")+geom_rug(color="gray50")

Again, you can see that the best-fitting line (dotted) has a slope somewhat less steep than before (compare to dashed line). Indeed, the slope is now r∑ZxZyn−1=3.55+3.898=7.448=0.93. Because the y-intercept is zero, the equation for the dotted line is Zy=.93×Zx. This will always be the case: Wen we have standardized variables, the slope of the line of best fit is equal to the correlation coefficient r.

OK, so what happens if our data are not linear? (We will go ahead and standardize our variables to avoid issues this time.)

data1<-data
data1$Y<--(data$Y^2)+1
data1<-data.frame(scale(data1))
ggplot(data1,aes(X,Y))+geom_point()+geom_hline(yintercept=0)+geom_vline(xintercept=0)+scale_x_continuous(limits=c(-3,3))+scale_y_continuous(limits = c(-3,3))+xlab("Zx")+ylab("Zy")+annotate("text",2,2,label=round(sum(data1[data1$X>0 & data1$Y>0,]$X*data1[data1$X>0 & data1$Y>0,]$Y),2))+
  annotate("text",-2,-2,label=round(sum(data1[data1$X<0 & data1$Y<0,]$X*data1[data1$X<0 & data1$Y<0,]$Y),2))+
  annotate("text",-2,2,label=round(sum(data1[data1$X<0 & data1$Y>0,]$X*data1[data1$X<0 & data1$Y>0,]$Y),2))+
  annotate("text",2,-2,label=round(sum(data1[data1$X>0 & data1$Y<0,]$X*data1[data1$X>0 & data1$Y<0,]$Y),2))+xlab("Zx")+ylab("Zy")+geom_rug(color="gray50")

#mean(data1$Y)
#sd(data1$Y)

We have r=∑ZxZyn−1=.65+(−.65)+2.61+(−2.61)8=08=0, zero correlation! This is because every point in a positive quadrant is cancelled out by a point in a negative quadrant!

However, there is still clearly a relationship between Zy and Zx! (Specifically, Zy=−Zx2+1). So calculating the correlation coefficient is only appropriate when you are interested in a linear relationship. As we will see, the same applies to linear regression.

Finally, let’s see what more realistic data might look like:

library(ggplot2)
library(MASS)

samples = 9
r = 0.5

data1 = as.data.frame(mvrnorm(n=samples, mu=c(0, 0), Sigma=matrix(c(1, r, r, 1), nrow=2), empirical=TRUE))
names(data1)<-c("X","Y")

ggplot(data1,aes(X,Y))+geom_point()+geom_hline(yintercept=0)+geom_vline(xintercept=0)+scale_x_continuous(limits=c(-3,3))+scale_y_continuous(limits = c(-3,3))+xlab("Zx")+ylab("Zy")+annotate("text",2,2,label=round(sum(data1[data1$X>0 & data1$Y>0,]$X*data1[data1$X>0 & data1$Y>0,]$Y),2))+
  annotate("text",-2,-2,label=round(sum(data1[data1$X<0 & data1$Y<0,]$X*data1[data1$X<0 & data1$Y<0,]$Y),2))+
  annotate("text",-2,2,label=round(sum(data1[data1$X<0 & data1$Y>0,]$X*data1[data1$X<0 & data1$Y>0,]$Y),2))+
  annotate("text",2,-2,label=round(sum(data1[data1$X>0 & data1$Y<0,]$X*data1[data1$X>0 & data1$Y<0,]$Y),2))+xlab("Zx")+ylab("Zy")+geom_rug(color="gray50")

Here, we find that r=∑ZxZyn−1=(−0.18)+3.1+(−1.01)+2.098=48=0.5

A positive correlation, but not nearly as strong as before.

To recap,

• r = strength of linear relationship between two variables

• for variables x and y, measures how close points (x,y) are to a straight line

• linear transformations (add/subtract a constant, multiply/divide by a constant) do not affect the value of r.

Now then, in the above case where the line fit perfectly (r=1), we could write our prediction for y (denoted ˆy) like so:

Zy=Zx⟶ˆy−ˉysdy=x−ˉxsdx

In this case, it was a perfect prediction! But when it is not perfect prediction, we need to calculate a slope coefficient for Zx

Regression

Correlation and regression are fundamental to modern statistics.

Regression techniques are useful when you want to model an outcome variable y as a function of (or in relationship to) predictor variables x1,x2,.... Specifically, regression is useful when you want

y=f(x1,x2,...)

Because this goal is so common, and because the tool is so general, regression models are ubiquitous in research today.

We are going to be talking here about situations where we model y as a linear function of x1,x2,...

y=b0+b1x1+b2x2+...+ϵ

We are interested in modeling the mean of our response variable y given some predictors x1,x2,.... Right now, think of regression as a function which gives us the expected value of y (ˆy) for a given value of x. The value ˆy is the predicted value, which will lie on the straight line that runs right through the middle of our scatterplot. Notice that some of the data will be above the predicted line and some will be below it: our ˆys will not perfectly match our y. The distances between the predicted (or “fitted”) values ˆy and the actual values y are called residuals (ϵ=y−ˆy). The residuals end up being very important conceptually: they represent the variability in y that is left after accounting for the best-fitting linear combination of your predictor(s) x (i.e., ˆy).

##Deriving regression coefficients!

Let’s keep it basic, because doing so will allow us to do something amazing using pretty basic math.

Let’s look at a special regression with one dependent variable y, one predictor x, and no intercept. Like we said above, the xs and ys we observed aren’t going to have a perfectly linear relationship: we need a residual term (ϵ) that capture deviations of our data from the predicted line ˆy=bx (“error”). Here’s the model:

y=bx+ϵˆy=bxsubtract bottom from top,y−ˆy=ϵ

Now, when we estimate b, our goal is to miminize the errors (the residuals ϵ=y−ˆy). Minimizing the errors/residuals will give us the best fitting line! But ϵi can be positive or negative, so we square them. Specifically, we want to choose b to minimize the sum of these squared residuals. Let’s do it! From above, we know that ϵ=y−bx, so ∑iϵ2i=∑i(yi−bxi)2. Let’s do a tiny bit of algebra (minb just means our goal is to find the value of b that minimizes the function:

minb∑iϵ2i=∑i(yi−ˆyi)2minb∑iϵ2i=∑i(yi−bxi)2=∑i(yi−bxi)(yi−bxi)=∑iy2i−2bxiyi+b2x2i

Recall that x and y are fixed. We want to know the value of b that will result in the smallest total squared error. We can find this value by taking the derivative of the error function (above) with respect to b: this results in a new function that tells you the slope of a line tangent to the original function for each value of b. Since the original function ∑i(yi−bxi)2 is a positive quadratic, we know that it has its minimum value at the bottom of the bend where the slope of the tangent line is zero. This is the minimum value of the error function! We just need to solve for b.

dropping subscripts for readabilityminb∑ϵ2=∑(y2−2bxy+b2x2)distribute sums, pull out constants, apply derivative operatorddb∑ϵ2=ddb∑y2−2b∑xy+b2∑x2=0−2∑xy+2b∑x2

set derivative equal to zero, solve for b−2∑xy+2b∑x2=02b∑x2=2∑xyb=∑xy∑x2

So the “least squares estimator” ˆb is ∑xy∑x2.

This means that our best prediction of ˆy given x is

ˆy=ˆbx=(∑ixiyi∑ix2i)x=∑iyixi∑ix2ix

So our prediction for y given x is a weighted average of all the observed yis.

This can be hard to see, so let’s say we observe the following data: y={−1,0,1},x={−1,1,0} Then given that x=1, our estimate ˆy is ∑iyixi∑ix2ix=−1−121+0121+1021=.5+0+0=.5

ggplot(data.frame(x=c(-1,1,0),y=c(-1,0,1)),aes(x,y))+geom_point()+geom_smooth(method=lm, se=F)+annotate(geom="point",x=1,y=.5,color="red")

Look again at ˆb=∑xy∑x2 for a moment. The top part of the estimator looks kind of like covariance, and the bottom looks kind of like a variance. Indeed,

cov(X,Y)var(X)=∑(xi−ˉx)(yi−ˉy)/n−1∑(xi−ˉx)2/n−1=∑(xi−ˉx)(yi−ˉy)∑(xi−ˉx)2

The only difference is that cov(X,Y)var(X) scales the variables first, by subtracting the mean from each! Remember that in our simplified equation yi=bxi+ϵi, we are implicitly assuming that the y-intercept is zero (as it is when both variables are standardized). When it is, you are fine to use ˆb=∑xiyi∑x2i However, when you do have a y-intercept you want to find (which is usually the case), then the best estimator of the true slope b1 is ˆb1=cov(X,Y)var(X)=∑(xi−ˉx)(yi−ˉy)∑(xi−ˉx)2. Let me show you the difference with some play data.

x1<-1:20
y1<-rnorm(20)+.5*x1+5
data<-data.frame(x1,y1)
b_noint=sum(x1*y1)/sum(x1^2)
b_int=cov(x1,y1)/var(x1)

ggplot(data,aes(x1,y1))+geom_point()+scale_x_continuous(limits=c(0,20))+scale_y_continuous(limits=c(0,15))+geom_abline(slope=b_noint,intercept=0)

y2<-scale(y1)
x2<-scale(x1)
ggplot(data,aes(x2,y2))+geom_point()+geom_abline(slope=b_noint,intercept=0)

I wanted to show a derivation in the simplest case, but it works just as well with more variables. Here is a run-through of the same, but with vectors and matricies. let y=Xb+e where y is a vector, X is a matrix of all predictors, b is a vector of parameters to be estimated, and e is a vector of errors.

[y1y2⋮]=[x1,1x1,2x2,1x2,2⋮⋮][b1b2⋮]+[ϵ1ϵ2⋮]

e′e=(y−Xb)′(y−Xb)e′e=y′y−2b′X′y+b′X′Xbddbe′e=−2Xy+2b′X′Xset equal to zero and solve for b0=−2X′y+2b′X′X2X′y=2b′X′Xb=(X′X)−1X′y

Even if you are relatively unfamiliar with matrix algebra, you can see the resemblance that (X′X)−1X′y bears to 1∑x2∑xy.

###Testing (X′X)−1X′Y on some data

Let’s say our data are Y=10,13,14 and X=2,5,12. Let’s derive the regression coefficient estimates ˆb=b0,b1 using the estimator X′X−1X′y

yxx=xxxxxxXbxxx+xxe[101314]=[1215112][b0b1]+[ϵ1ϵ2ϵ3]

Now ˆb=(X′X)−1X′Y, so... Click button to expand the calculation: Show/Hide

So we found that the intercept estimate is 10.09 and the slope estimate is 0.35. Let’s check this result by conducting a regression on the same data:

lm(c(10,13,14)~c(2,5,12))

## 
## Call:
## lm(formula = c(10, 13, 14) ~ c(2, 5, 12))
## 
## Coefficients:
## (Intercept)  c(2, 5, 12)  
##     10.0886       0.3544

dat<-data.frame(x=c(2,5,12),y=c(10,13,14))

lm_eqn <- function(df){
    m <- lm(y ~ x, df);
    eq <- substitute(hat(italic(y)) == a + b %.% italic(x), 
         list(a = format(coef(m)[1], digits = 2), 
              b = format(coef(m)[2], digits = 2)))
    as.character(as.expression(eq));                 
}

ggplot(dat,aes(x,y))+geom_point()+geom_smooth(method=lm,se=F)+geom_text(x = 5, y = 14, label = lm_eqn(dat), parse = TRUE)