Return to Answer

deleted 138 characters in body

Source Link

edited Feb 10, 2023 at 17:29

82.8k
32
201
663

x <- rep(1:10, each=2)
y <- x/10 + 0:1
plot(x,y)

sad <- function(x,y,coef) { # mad is sad/n
    yhat <- coef[1] + coef[2]*x
    resid <- y - yhat
    sum( abs( resid ) )
}

library(quantreg)
fit0 <- rq( y~x )
abline(fit0)

fit1 <- lm( y~x, subset= c(1,20) )
fit2 <- lm( y~x, subset= c(2,19) )
fit3 <- lm( y~x, subset= c(2,20) )
fit4 <- lm( y~x, subset= c(1,19) )

fit5.coef <- c(0.5, 1/10)

abline(fit1)
abline(fit2)
abline(fit3)
abline(fit4)
abline(fit5.coef)
for (i in seq( -0.5, 0.5, by=0.1 ) ) {
    abline( fit5.coef + c(i,0) )
}

tmp1 <- seq( coef(fit1)[1], coef(fit2)[1], len=10 )
tmp2 <- seq( coef(fit1)[2], coef(fit2)[2], len=10 )

for (i in seq_along(tmp1) ) {
    abline( tmp1[i], tmp2[i] )
}

sad(x,y, coef(fit0))
sad(x,y, coef(fit1))
sad(x,y, coef(fit2))
sad(x,y, coef(fit3))
sad(x,y, coef(fit4))
sad(x,y, fit5.coef )

for (i in seq( -0.5, 0.5, by=0.1 ) ) {
    print(sad(x,y, fit5.coef + c(i,0) ))
}

for (i in seq_along(tmp1) ) {
    print(sad(x,y, c(tmp1[i], tmp2[i]) ) )
}

set.seed(1)
y2 <- y + rnorm(20,0,0.25)
plot(x,y2)
fitnew <- rq(y2~x)  # note the still non-unique warning
abline(fitnew)
abline(coef(fitnew) + c(.1,0))
abline(coef(fitnew) + c(0, 0.01) )
sad( x,y2, coef(fitnew) )
sad( x,y2, coef(fitnew)+c(.1,0))
sad( x,y2, coef(fitnew)+c(0,0.01))

x <- rep(1:10, each=2)
y <- x/10 + 0:1
plot(x,y)
    
sad <- function(x,y,coef) { # mad is sad/n
        yhat <- coef[1] + coef[2]*x
        resid <- y - yhat
        sum( abs( resid ) )
    }
    
library(quantreg)
fit0 <- rq( y~x )
abline(fit0)
    
fit1 <- lm( y~x, subset= c(1,20) )
fit2 <- lm( y~x, subset= c(2,19) )
fit3 <- lm( y~x, subset= c(2,20) )
fit4 <- lm( y~x, subset= c(1,19) )
    
fit5.coef <- c(0.5, 1/10)
    
abline(fit1)
abline(fit2)
abline(fit3)
abline(fit4)
abline(fit5.coef)
for (i in seq( -0.5, 0.5, by=0.1 ) ) {
        abline( fit5.coef + c(i,0) )
    }
    
tmp1 <- seq( coef(fit1)[1], coef(fit2)[1], len=10 )
tmp2 <- seq( coef(fit1)[2], coef(fit2)[2], len=10 )
    
for (i in seq_along(tmp1) ) {
        abline( tmp1[i], tmp2[i] )
    }
    
sad(x,y, coef(fit0))
sad(x,y, coef(fit1))
sad(x,y, coef(fit2))
sad(x,y, coef(fit3))
sad(x,y, coef(fit4))
sad(x,y, fit5.coef )
    
for (i in seq( -0.5, 0.5, by=0.1 ) ) {
        print(sad(x,y, fit5.coef + c(i,0) ))
    }
    
for (i in seq_along(tmp1) ) {
        print(sad(x,y, c(tmp1[i], tmp2[i]) ) )
    }
    
set.seed(1)
y2 <- y + rnorm(20,0,0.25)
plot(x,y2)
fitnew <- rq(y2~x)  # note the still non-unique warning
abline(fitnew)
abline(coef(fitnew) + c(.1,0))
abline(coef(fitnew) + c(0, 0.01) )
sad( x,y2, coef(fitnew) )
sad( x,y2, coef(fitnew)+c(.1,0))
sad( x,y2, coef(fitnew)+c(0,0.01))

x <- rep(1:10, each=2)
y <- x/10 + 0:1
plot(x,y)

sad <- function(x,y,coef) { # mad is sad/n
    yhat <- coef[1] + coef[2]*x
    resid <- y - yhat
    sum( abs( resid ) )
}

library(quantreg)
fit0 <- rq( y~x )
abline(fit0)

fit1 <- lm( y~x, subset= c(1,20) )
fit2 <- lm( y~x, subset= c(2,19) )
fit3 <- lm( y~x, subset= c(2,20) )
fit4 <- lm( y~x, subset= c(1,19) )

fit5.coef <- c(0.5, 1/10)

abline(fit1)
abline(fit2)
abline(fit3)
abline(fit4)
abline(fit5.coef)
for (i in seq( -0.5, 0.5, by=0.1 ) ) {
    abline( fit5.coef + c(i,0) )
}

tmp1 <- seq( coef(fit1)[1], coef(fit2)[1], len=10 )
tmp2 <- seq( coef(fit1)[2], coef(fit2)[2], len=10 )

for (i in seq_along(tmp1) ) {
    abline( tmp1[i], tmp2[i] )
}

sad(x,y, coef(fit0))
sad(x,y, coef(fit1))
sad(x,y, coef(fit2))
sad(x,y, coef(fit3))
sad(x,y, coef(fit4))
sad(x,y, fit5.coef )

for (i in seq( -0.5, 0.5, by=0.1 ) ) {
    print(sad(x,y, fit5.coef + c(i,0) ))
}

for (i in seq_along(tmp1) ) {
    print(sad(x,y, c(tmp1[i], tmp2[i]) ) )
}

set.seed(1)
y2 <- y + rnorm(20,0,0.25)
plot(x,y2)
fitnew <- rq(y2~x)  # note the still non-unique warning
abline(fitnew)
abline(coef(fitnew) + c(.1,0))
abline(coef(fitnew) + c(0, 0.01) )
sad( x,y2, coef(fitnew) )
sad( x,y2, coef(fitnew)+c(.1,0))
sad( x,y2, coef(fitnew)+c(0,0.01))

x <- rep(1:10, each=2)
y <- x/10 + 0:1
plot(x,y)
    
sad <- function(x,y,coef) { # mad is sad/n
        yhat <- coef[1] + coef[2]*x
        resid <- y - yhat
        sum( abs( resid ) )
    }
    
library(quantreg)
fit0 <- rq( y~x )
abline(fit0)
    
fit1 <- lm( y~x, subset= c(1,20) )
fit2 <- lm( y~x, subset= c(2,19) )
fit3 <- lm( y~x, subset= c(2,20) )
fit4 <- lm( y~x, subset= c(1,19) )
    
fit5.coef <- c(0.5, 1/10)
    
abline(fit1)
abline(fit2)
abline(fit3)
abline(fit4)
abline(fit5.coef)
for (i in seq( -0.5, 0.5, by=0.1 ) ) {
        abline( fit5.coef + c(i,0) )
    }
    
tmp1 <- seq( coef(fit1)[1], coef(fit2)[1], len=10 )
tmp2 <- seq( coef(fit1)[2], coef(fit2)[2], len=10 )
    
for (i in seq_along(tmp1) ) {
        abline( tmp1[i], tmp2[i] )
    }
    
sad(x,y, coef(fit0))
sad(x,y, coef(fit1))
sad(x,y, coef(fit2))
sad(x,y, coef(fit3))
sad(x,y, coef(fit4))
sad(x,y, fit5.coef )
    
for (i in seq( -0.5, 0.5, by=0.1 ) ) {
        print(sad(x,y, fit5.coef + c(i,0) ))
    }
    
for (i in seq_along(tmp1) ) {
        print(sad(x,y, c(tmp1[i], tmp2[i]) ) )
    }
    
set.seed(1)
y2 <- y + rnorm(20,0,0.25)
plot(x,y2)
fitnew <- rq(y2~x)  # note the still non-unique warning
abline(fitnew)
abline(coef(fitnew) + c(.1,0))
abline(coef(fitnew) + c(0, 0.01) )
sad( x,y2, coef(fitnew) )
sad( x,y2, coef(fitnew)+c(.1,0))
sad( x,y2, coef(fitnew)+c(0,0.01))

Source Link

created Jan 30, 2013 at 19:11

Greg Snow

52.9k
2
112
190

This is more a response to @PeterFlom's comment on my comment, but it is too big to fit in a comment (and does relate to the original question).

Here is some R code to show a case where there are multiple lines that all give the same minimum MAD/SAD values.

The first part of the example is clearly contrived data to demonstrate, but the end includes more of a random element to demonstrate that the general concept will still hold in some more realistic cases.

x <- rep(1:10, each=2)
y <- x/10 + 0:1
plot(x,y)

sad <- function(x,y,coef) { # mad is sad/n
    yhat <- coef[1] + coef[2]*x
    resid <- y - yhat
    sum( abs( resid ) )
}

library(quantreg)
fit0 <- rq( y~x )
abline(fit0)

fit1 <- lm( y~x, subset= c(1,20) )
fit2 <- lm( y~x, subset= c(2,19) )
fit3 <- lm( y~x, subset= c(2,20) )
fit4 <- lm( y~x, subset= c(1,19) )

fit5.coef <- c(0.5, 1/10)

abline(fit1)
abline(fit2)
abline(fit3)
abline(fit4)
abline(fit5.coef)
for (i in seq( -0.5, 0.5, by=0.1 ) ) {
    abline( fit5.coef + c(i,0) )
}

tmp1 <- seq( coef(fit1)[1], coef(fit2)[1], len=10 )
tmp2 <- seq( coef(fit1)[2], coef(fit2)[2], len=10 )

for (i in seq_along(tmp1) ) {
    abline( tmp1[i], tmp2[i] )
}

sad(x,y, coef(fit0))
sad(x,y, coef(fit1))
sad(x,y, coef(fit2))
sad(x,y, coef(fit3))
sad(x,y, coef(fit4))
sad(x,y, fit5.coef )

for (i in seq( -0.5, 0.5, by=0.1 ) ) {
    print(sad(x,y, fit5.coef + c(i,0) ))
}

for (i in seq_along(tmp1) ) {
    print(sad(x,y, c(tmp1[i], tmp2[i]) ) )
}

set.seed(1)
y2 <- y + rnorm(20,0,0.25)
plot(x,y2)
fitnew <- rq(y2~x)  # note the still non-unique warning
abline(fitnew)
abline(coef(fitnew) + c(.1,0))
abline(coef(fitnew) + c(0, 0.01) )
sad( x,y2, coef(fitnew) )
sad( x,y2, coef(fitnew)+c(.1,0))
sad( x,y2, coef(fitnew)+c(0,0.01))