On Feb 11, 2019, at 10:24 AM, Marco Besozzi <marco.beso48 at gmail.com> wrote:
I employed the "galton" set of data included in the package "psych". With
the package "mblm" I obtained the Theil-Sen nonparametric regression and
the Siegel non parametric regression, and compared them with the ordinary
least square regression line.
The results of standard regression and Theil-Sen regression are practically
identical. But the Siegel regression seems to have a bias that I cannot
understand. May I ask for a possible explanation? The bias may be related
to the number of ties in the set of data? Here's the code and the image.
Best regards.
Marco Besozzi
# Theil-Sen and Siegel nonparametric regression with package mblm
# comparison with ordinary least squares (parametric) regression
# on galton set of data included in the package psych
#
library(psych)
attach(galton)
library(mblm)
#
reglin_yx <- lm(child ~ parent, data=galton) # ordinary least squares
(parametric) regression
a_yx <- reglin_yx$coefficients[1] # intercept a
b_yx <- reglin_yx$coefficients[2] # slope b
#
regnonTS <- mblm(child ~ parent, data=galton, repeated=FALSE) # Theil-Sen
nonparametric regression (wait a few minutes!)
a_TS <- regnonTS$coefficients[1] # intercept a
b_TS <- regnonTS$coefficients[2] # slope b
#
regnonS = mblm(child ~ parent, data=galton, repeated=TRUE) # Siegel
nonparametric regression
a_S <- regnonS$coefficients[1] # intercept a
b_S <- regnonS$coefficients[2] # slope b
#
# xy plot of data and regression lines
#
windows() # open a new window
plot(parent, child, xlim = c(60,80), ylim = c(60,80), pch=1, xlab="Parent
heigt (inch)", ylab="Chile height (inch)", main="Regression lines
comparison", cex.main = 0.9) # data plot
abline(a_yx, b_yx, col="green", lty=1) # ordinary least squares
(parametric) regression line
abline(a_TS, b_TS, col="blue", lty=1) # Theil-Sen nonparametric regression
line
abline(a_S, b_S, col="red", lty=1) # Siegel nonparametric regression
legend(60, 80, legend=c("Ordinary least squares regression", "Theil-Sen
nonparametric regression","Siegel nonparametric regression"),
col=c("green", "blue", "red"), lty=c(4,4,1), cex=0.8) # add a legend
#
<Siegel.PNG>______________________________________________