Return to Answer

#set up some example data
alpha<-1.4
beta<-0.87
x<-runif(20,3,15)
epsilon<-rlogis(20,0,0.8)
y<-alpha+beta*x+epsilon

# start with a robust line (can substitute lm here if you want)
#   to get a rough idea of an interval. This can be improved!
ab<-line(x,y)$coefficients
bi<-ab[2]

# The next two lines of code are the actual algorithm:
#  (you may want to play with some of the other arguments to uniroot)
T<-function(b,y,x) {xm <- mean(x); sum((x-xm)*rank(y-b*x))}
Tzero<-uniroot(T,interval=sort(c(0,2*bi)),y=y,x=x,extendInt="yes")

b<-Tzero$root   #pull out the coefficient

#---- illustration of performance in four plots
op<-par()
par(mfrow=c(2,2))

#1 don't use this normally, just to illustrate that it works
bb<-seq(0,2*b,l=101)
tt<-array(NA,length(bb))
for (i in seq_along(bb)) tt[i]=T(bb[i],y,x)  
plot(tt~bb,pch=16,cex=0.1,xlab="beta",ylab="T",main="Hodges-Lehmann T function" )
abline(h=0,v=b,col=8,lty=3)

#2 
plot(x,y,main="Data with fitted line")
a<-median(y-b*x) # substitute whatever intercept you like
abline(a=a,b=b,col=2)

#3
fitted<-a+b*x
residual<-y-fitted
plot(fitted,residual,main="residuals vs fitted values")
abline(h=0,col=8)

#4
hist(residual)

#restore plot parameters
par(op)

    #set up some example data
    alpha<-1.4
    beta<-0.87
    x<-runif(20,3,15)
    epsilon<-rlogis(20,0,0.8)
    y<-alpha+beta*x+epsilon

    # start with a robust line (can substitute `lm` here if you want)
    #   to get a rough idea of an interval. This can be improved!
    ab<-line(x,y)$coefficients
    bi<-ab[2]

    # The next two lines of code are the actual algorithm:
    #  (you may want to play with some of the other arguments 
    #  to uniroot)
    T <- function(b,y,x) {xm <- mean(x); sum((x-xm)*rank(y-b*x))}
    Tzero <- uniroot(T, interval=sort(c(0,2*bi)), y=y, x=x, 
                   extendInt="yes")
    
    b <- Tzero$root   #pull out the coefficient

    #---- illustration of performance in four plots
    op<-par()
    par(mfrow=c(2,2))
    
    #1 don't use this normally, just to illustrate that it works
    bb<-seq(0,2*b,l=101)
    tt<-array(NA,length(bb))
    for (i in seq_along(bb)) tt[i]=T(bb[i],y,x)  
    plot(tt~bb, pch=16, cex=0.1, xlab="beta", ylab="T", 
         main="Hodges-Lehmann T function" )
    abline(h=0, v=b, col=8, lty=3)
    
    #2 
    plot(x,y,main="Data with fitted line")
    a<-median(y-b*x) # substitute whatever intercept you like
    abline(a=a,b=b,col=2)
    
    #3
    fitted<-a+b*x
    residual<-y-fitted
    plot(fitted,residual,main="residuals vs fitted values")
    abline(h=0,col=8)
    
    #4
    hist(residual)
    
    #restore plot parameters
    par(op)

A suitable algorithm will require you to find a good starting point, and many will also require you to start by bracketing the root, which may require an initial search procedure. I sort of dodge the search for a bracket in the example below because the function I call will extend the bracket if it's insufficient (we could improve that part substantially).

Note that the root finder stops with some tolerance on the root (which you have control over). Because the function has jumps, the value of the function at the returned root may not seem all that "close" to 0, but it will jump to the other side of 0 with only a small shift.

# The next two lines of code are the actual algorithm:
#  (you may want to play with some of the arguments to uniroot)
T<-function(b,y,x) {xm <- mean(x); sum((x-xm)*rank(y-b*x))}
Tzero<-uniroot(T,interval=sort(c(0,2*bi)),y=y,x=x,extendInt="yes")

b<-Tzero$root   #pull out the coefficient

Note that the slope for the "Spearman line" in the answer here was identified by root-finding (as was the quadrant-correlation slope mentioned near the end).

A suitable algorithm will require you to find a good starting point, and/or will require you to start by bracketing the root, which bracketing may require an initial search procedure. I sort of dodge the search for a bracket in the example below because the function I call will extend the bracket if it's insufficient (we could improve the bracketing part substantially).

Note that the root finder stops with some tolerance on the root (which you have control over in the root finding function). Because the function has jumps, the value of the function at the returned root may not seem all that "close" to 0, but it will jump to the other side of 0 with only a small shift.

# The next two lines of code are the actual algorithm:
#  (you may want to play with some of the other arguments to uniroot)
T<-function(b,y,x) {xm <- mean(x); sum((x-xm)*rank(y-b*x))}
Tzero<-uniroot(T,interval=sort(c(0,2*bi)),y=y,x=x,extendInt="yes")

b<-Tzero$root   #pull out the coefficient

Note that the slope for the "Spearman line" in the answer here was also identified by root-finding (as was the quadrant-correlation slope mentioned near the end).

I haven't written this as a full function of its own, but it's easy enough to wrap in a function.

# The next two lines of code are the actual algorithm:
#  (you may want to play with some of the arguments to uniroot)
T<-function(b,y,x) {xm <- mean(x); sum((x-xm)*rank(y-b*x))}
Tzero<-uniroot(T,sort(c(0,2*bi)),y=y,x=x,extendInt="yes")

b<-Tzero$root   #pull out the coefficient

I haven't written this as a full function of its own - the point was just to illustrate the basic idea was implemented as I said above - but it's easy enough to wrap in a function.

# The next two lines of code are the actual algorithm:
#  (you may want to play with some of the arguments to uniroot)
T<-function(b,y,x) {xm <- mean(x); sum((x-xm)*rank(y-b*x))}
Tzero<-uniroot(T,interval=sort(c(0,2*bi)),y=y,x=x,extendInt="yes")

b<-Tzero$root   #pull out the coefficient