# R: Regression, Two-stage Least Squares

I need to perform manually two-stage Least Squares(to illustrate its advantages), where the first stage is repeated median estimate and the second stage should be weighted least squares, where weights are obtained(as far, as I understand) from polynomial regression of first-stage residuals on regressors.

Suppose I have generated the following heteroscedastic model:

where error depends on regressor:

set.seed(100)
b<-c(12,7.25) ## my coefficients
num<-50 ## number of observations

raw_x<-runif(num,min=0,max=2) ## regressors

my_y<-as.vector(b%*%t(data.frame(rep(1,num),raw_x))+
rnorm(num,mean=0,sd=(raw_x-1)^2)) ## observations

l<-lm(my_y~raw_x) ## let's create linear model

plot(fitted(l),residuals(l)) ## we see heteroskedasticity
## we got to higher values, our residuals explode

abline(0,0)
title("Residual vs Fit. value");


So I perform repeated median regression (formula in the Introduction):

## Generating first model using repeated median

## slope

fij = function(i,j)
{
(my_y[i]-my_y[j])/(raw_x[i]-raw_x[j])
}

bij<-outer(1:num,1:num,fij) ##NaN's were produced on the diagonal

rowmeds <- apply(bij, 1, median,na.rm=TRUE)
b_med<-median(rowmeds)

# colmeds <- apply(bij, 2, median,na.rm=TRUE) ## column medians are the same
# b_med3<-median(colmeds)

## Intercept

## med(y_i - b*x_i)

a_med<-median(my_y-b_med*raw_x)


The fit is extremely accurate! In this example a_med is 11.97634 and b_med equals 7.27022.

Now I perform 2nd order polynomial regression of residuals:

so that (X here is m x 3 matrix):

I was told that as long as residual variances can be roughly estimated from only one observation, actual fit from this model can be used; residual variances = coefficients for the weighted least squares:

## Obtaining 2nd order polynomial estimator for residuals

Xmatr = t(rbind(rep(1,num),raw_x,(raw_x)^2))

## coef_var = (X^T*X)^(-1)*X^T^e
coef_var<-solve(t(Xmatr)%*%Xmatr)%*%t(Xmatr)%*%t(my_y-c(a_med,b_med)%*%rbind(rep(1,num),raw_x))

## Obtaining sigma (residual deviation) esitmate

my_sigma<-t(Xmatr%*%coef_var)

## performing regression with sigma-weights

b_wls<-lm(as.numeric(my_y)~raw_x,weights=as.numeric(my_sigma^2))$coef ## final plot library(scales) plot(raw_x,my_y, pch=20,col=alpha("salmon",0.6)) abline(b[1],b[2], col="black") ## real line abline(a_med,b_med,col="blue") ## repeated median fit abline(b_wls, col="magenta") legend('bottomright', c("Real","Repeat median","Two-Level LS") , lty=1, col=c('black', 'blue','magenta'), bty='n', cex=.75)  The resulting fit is always worse (and sometimes turns into complete garbage). Please, can you explain me what I'm doing wrong? I need to obtain the result where two-level LS is better than repeated median fit(provided the error depends on regressors as shown before). EDIT: Using R function lm seems to produce the same picture: Xmatr = t(rbind(rep(1,num),raw_x,(raw_x)^2)) residual<-as.vector(my_y-c(a_med,b_med)%*%rbind(rep(1,num),raw_x)) polycoef<-lm(residual ~ poly(raw_x, 2, raw=TRUE))$coefficients
my_sigma<-t(Xmatr%*%polycoef)


EDIT2: Checking the quality of the residual fit (as far as I understand what's going on):

plot(raw_x,abs(residual))
lines(sort(raw_x),(sort(raw_x)-1)^2) ## real residuals
lines(sort(raw_x), my_sigma[order(raw_x)],col = "magenta") ## fitted residuals
legend('topleft', c("Real res","Fitted res") ,
lty=1, col=c('black','magenta'), bty='n', cex=.75)


EDIT3: As long as my standard deviation is $(X_i-1)^2$, so the variance has power 4 and maybe I should do for the residual fit

residual<-abs(as.vector(my_y-c(a_med,b_med)%*%rbind(rep(1,num),raw_x)))


this makes residual fit by the quadratic function better, but the regression line moves even more far than before.

## migrated from stackoverflow.comOct 20 '15 at 14:26

This question came from our site for professional and enthusiast programmers.

Particular solution: looks like "weights" parameter in lm() is the diagonal of $W^{-1}$ in the formula of weighted LS: $b_{wls} = (X^T W^{-1} X)^{-1}W^{-1}y$, $W = diag(\sigma_1,...\sigma_m)$. So
b_wls<-lm(as.numeric(my_y)~raw_x,weights=as.numeric(1/my_sigma^2))\$coef

raw_x<-runif(num,min=0,max=9)