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$

closed as off-topic by mkt - Reinstate Monica, Michael Chernick, Peter Flom - Reinstate Monica Oct 1 at 10:48

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – mkt - Reinstate Monica, Michael Chernick, Peter Flom - Reinstate Monica
If this question can be reworded to fit the rules in the help center, please edit the question.

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$

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