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
#
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Siegel.PNG
Type: image/png
Size: 81522 bytes
Desc: not available
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20190211/ac668c7b/attachment.png>