4
$\begingroup$

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:

Y_i = b_0+b_1 X_i + \epsilon_i

where error depends on regressor:

\epsilon_i \sim N(0,(X_i-1)^2)

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");

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:

\hat{\epsilon_i}=a_0 + a_1X_i+a_2X_i^2+\delta_i

\hat{\epsilon} = \begin{pmatrix} 1 & X_1 & X_1^2 \\ \vdots \\ 1 & X_m & X_m^2 \end{pmatrix}\begin{pmatrix}a_0 \\ a_1 \\ a_2 \end{pmatrix}+\delta

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

\hat{a} = [X^TX]^{-1}X^T\hat{\epsilon}

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:

\hat{\hat{\sigma}}=X\hat{a}

## 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).

enter image description here enter image description here

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)

enter image description here

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. enter image description here

$\endgroup$

migrated from stackoverflow.com Oct 20 '15 at 14:26

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

0
$\begingroup$

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

greatly improved the accuracy, but it is still not as good as repeated median in general case. If max value is big in

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

The picture crashes again.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.