Calculating Bare Soil Line using R?

Question

I am using a package in R called BSL (Bare soil line), package is at https://rdrr.io/cran/landsat/man/BSL.html. I've used it successfully to run regression line for my soil study area but it is Model II regression (specifically the lmodel2 package in R) and I don't understand why it's Model II regression and not a simple linear regression using the lm() function in R. Can someone explain why it would be using Model II regression instead of Simple Linear Regression?

A bit more info: the soil line is a linear relationship between reflectance values between Red and Near-infrared wavelengths and is defined by the equation:

NIR = αRed + β where NIR and Red correspond to the near-infrared and red bands of the satellite sensor, alpha is the slope, and beta the y intercept.

For my research, I'm discovering that more fertile soil may have lower y-intercept value and smaller slope values, but this seems to be based on the Model 2 regression values.

For my soil study area, I run the BSL command in R and I get this output:

> result.bsl$BSL
  Intercept       Slope
1323.007184    0.640505
> result.bsl$summary

Model II regression

Call: lmodel2(formula = bsl.joint[ratio43 < quantile(ratio43, llimit), 2]~ bsl.joint[ratio43 < quantile(ratio43, llimit), 1])

n = 3   r = 0.9878535   r-square = 0.9758546
Parametric P-values:   2-tailed = 0.09932547    1-tailed = 0.04966273
Angle between the two OLS regression lines = 0.6373716 degrees

Regression results
  Method Intercept     Slope Angle (degrees) P-perm (1-tailed)
1    OLS  1350.359 0.6359507        32.45435                NA
2     MA  1323.007 0.6405050        32.63976                NA
3    SMA  1303.397 0.6437702        32.77223                NA

Confidence intervals
  Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
1    OLS      -6286.498        8987.216 -0.6351029    1.907004
2     MA      48805.536        9004.094 -0.6384682   -7.265783
3    SMA     -11020.651        4246.396  0.1537332    2.695840

Eigenvalues: 66506.64 336.0233

H statistic used for computing C.I. of MA: 0.8240161

If I run lm(nir~red, mydata), I get a different slope and intercepts for the same data:

Call:
lm(formula = mydata$X1.5 ~ mydata$X1.3, data = mydata)

Residuals:
    Min      1Q  Median      3Q     Max 
-713.28 -100.97    9.46  118.72  505.87 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 454.07689  240.22340    1.89   0.0594  
mydata$X1.3   0.89429    0.03846   23.25   <2e-16 
Signif. codes:  0     0.001   0.01   0.05 . 0.1  1

Residual standard error: 185.8 on 439 degrees of freedom
Multiple R-squared:  0.5518,    Adjusted R-squared:  0.5508 
F-statistic: 540.6 on 1 and 439 DF,  p-value: < 2.2e-16

In the model 2 regression, my y intercept and slope are 1323 and 0.64 respectively. But in the lm model, my y intercept and slope are 454.08 and 0.89 respectively. I've studied statistics at the first year college level, so I seek any explanations about why model 2 Regression is used here.

Answer 1

Here's your model:

NIR = αRed + β where NIR and Red correspond to the near-infrared and red bands of the satellite sensor, alpha is the slope, and beta the y intercept.

An ordinary least squares model is that your explanatory variable (here, Red) is known with 0 error, and your response (NIR) has independent Normal errors with the same variance.

In R you'd fit this with lm(NIR ~ Red) and get slope and intercept. You might think that if you turned this around and fitted lm(Red ~ NIR) you'd get the same model fit but with everything reflected in the X-Y diagonal. But because the mathematics has that the explanatory variable has no error in it, you don't get that. The model is not symmetric under exchange of NIR and Red because ordinary least squares was developed for experiments where you can set one thing precisely (your explanatory variable) and measure another thing (the response) with error.

For the two satellite bands its not obvious that one is an "explanatory variable" and the other a "response". Both variables are equivalent, and both have sources of uncertainty. In this case the models fitted by lmodel2 consider both variables as having uncertainty.

It follows then that if you use lmodel2 and flip the formula round, you should get a model that is essentially the same model reflected in the X=Y diagonal. Let's try that with the sample model in the documentation:

> lmodel2(Survival ~ Predicted_by_model , data=mod2ex1, nperm=99)$regression.results
RMA was not requested: it will not be computed.

  Method   Intercept    Slope Angle (degrees) P-perm (1-tailed)
1    OLS -0.09110673 1.069598        46.92606              0.01
2     MA -0.65028338 1.334739        53.15907              0.01
3    SMA -0.52483724 1.275257        51.89803                NA
> lmodel2(Predicted_by_model ~ Survival , data=mod2ex1, nperm=99)$regression.results
RMA was not requested: it will not be computed.

  Method Intercept     Slope Angle (degrees) P-perm (1-tailed)
1    OLS 0.6852956 0.6576961        33.33276              0.01
2     MA 0.4871990 0.7492103        36.84093              0.01
3    SMA 0.4115541 0.7841557        38.10197                NA

For the OLS model the slopes are 46 and 33 degrees, which adds up to 79. For the two models specific to the lmodel2 package, the slope angles add up to 90 degrees, showing they are reflected in the X-Y diagonal.