Your LAD function:

LADf=function(par) {(sum(y-par[1]-par[2]*x)^2)}

Looks exactly like least squares for me. Hence what you do is minimizing sum of squares, not abs of deviance. You need something like

LADf <- function(par) { sum(abs(y - par[1] - par[2]*x)) }

Note that this function is not differentiable, so you have to use an optimizer that can handle that (Nelder-Mead or SANN for instance).

LAD estimator is also equivalent to median regression so you can do this with quantile regression package instead.

The definition of "analogue" is not clear.

If you view $R^2$ as a metric to measure the "goodness of the fit" in regression setting. Then, likelihood can be used to evaluate the "goodness of the fit" for LDA (Linear Discriminant Analysis).

EDIT: I think OP was trying to ask LAD (Least absolute deviations regression) but not LDA. Here is my answer to LAD.

I would suggest OP to review the loss function in regression setting. For squared loss regression, the loss function is

$$\sum_i (y_i-\hat y_i)^2$$

where $\hat y_i$ is the predicted value for data $i$, and the sum is over all the data points.

For least absolute deviation loss, the loss function is

$$\sum_i |(y_i-\hat y_i)|$$

Therefore, to evaluate the "goodness of fit", we can exam the loss value. Now, the problem is with the loss value metric, it will depend on number of data points, i.e., the more data points we have, in general, the higher value the will be in aforementioned 2 formulas.

One way of dealing it is divided by number of data points. And if we want to do one step more, we can "normalize" it into 0 to 1 using RSS, which is R square.

But I think dividing the number of data points should be good for your purpose.

Someone's going to complain that this question belongs to other stackoverflow forum. But the quick take on your question is, you need to follow the standard R idiom. That is, for function like lad, it wants a formula like y ~ x1 + ... and a data object. To use your syntax, you want l1fit which accepts a x and y. Just read the documentation carefully.

Your attempt didn't work because you specified "qr" as the method for geom_smooth rather than "rq", which gives you a least absolute deviation best fit line.

ggplot(data, aes(x=country_mean_rep, y=country_mean_crime)) + 
  geom_point() + 
  geom_smooth(aes(colour="Linear", fill="Linear"), 
              method="lm", 
              formula=y ~ x,
              se = F) +
  geom_smooth(aes(colour="Least Absolute Deviation", fill="Least Absolute Deviation"), 
              method="rq",
              se = F,
              formula = y ~ x) +
  labs(colour="Functional Form", fill="Functional Form") +
  geom_text(aes(label=country_name), nudge_y=0.02) +
  theme_bw()

Output:

The caveat of this solution is that for it to work, you have to set se = F for your least absolute deviation line, which means you can't plot confidence intervals for this model. You can still have your confidence interval for the linear model if you like; just set se = T in the respective geom_smooth call.

There's not actually a t-test, because the estimate divided by its standard error doesn't have a t-distribution. Similar for an F-test.

Being a maximum likelihood estimate, there would be an asymptotic z-test, or an asymptotic chi-square test.

[There's the possibility of using some resampling-based tests as well, permutation tests or bootstrapping. You could also use L1pack's ability to simulate from L1 models to do a parametric bootstrap.]

You could compute something like an $R^2$ but it doesn't quite make sense because you're not computing something that tries to maximize that; it might make more sense to compute some analogous statistic, but $R^2$ has a number of properties and it depends on which things you try to carry over and which you don't.

[I note that quantreg::rq (which by default does L1 regression) will give an interval for the coefficients; this allows for a test (since you can see if the interval includes 0). There are also some other testing options in that package]