Sometimes you have two variables that are correlated and where you can use x to predict y or y to predict x but there is no clear pattern of causality between the two. This might be two biological measurements that share a common cause or ir might be two methods for measuring the same value.
Problems like this are common in allometry and in methods research. If you want to carry out a regression in these cases you should be using a Deming regression. You can do this in R with the Deming library.
I am going to use a subset of the Iris data for a single species as an example analysis.
I can plot the OLS and Standardised Major Axis model for the setosa species data.
library(deming)
library(ggplot2)
library(lmodel2)
library(ggpmisc)
data(iris)
setosa <- subset(iris, Species=="setosa")
ggplot(subset(setosa), aes(x=Petal.Length, y=Petal.Width))+
geom_point()+
geom_smooth(method=lm, colour="#4b0082")+
stat_poly_eq(mapping = use_label("eq", "R2", "n", "P"),colour="#4b0082")+
stat_ma_line(method="SMA",colour="darkorange") +
stat_ma_eq(mapping = use_label("eq", "R2", "n", "P"), orientation = "y",label.y = 0.9, colour="darkorange")
First I have calculated the OLS and Standardised Major Axis models. Then I can compare this to the Deming model.
library(deming)
fit1.setosa <- lm(Petal.Width~Petal.Length, setosa)
summary(fit1.setosa)
Call:
lm(formula = Petal.Width ~ Petal.Length, data = setosa)
Residuals:
Min 1Q Median 3Q Max
-0.15365 -0.05365 -0.03352 0.06632 0.32623
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.04822 0.12164 -0.396 0.6936
Petal.Length 0.20125 0.08263 2.435 0.0186 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1005 on 48 degrees of freedom
Multiple R-squared: 0.11, Adjusted R-squared: 0.09144
F-statistic: 5.931 on 1 and 48 DF, p-value: 0.01864confint(fit1.setosa) 2.5 % 97.5 %
(Intercept) -0.29279626 0.1963556
Petal.Length 0.03510125 0.3673889fit2.setosa <- lmodel2(Petal.Width~Petal.Length, setosa, nperm=999)
fit2.setosa
Model II regression
Call: lmodel2(formula = Petal.Width ~ Petal.Length, data = setosa,
nperm = 999)
n = 50 r = 0.33163 r-square = 0.1099785
Parametric P-values: 2-tailed = 0.01863892 1-tailed = 0.009319458
Angle between the two OLS regression lines = 49.96531 degrees
Permutation tests of OLS, MA, RMA slopes: 1-tailed, tail corresponding to sign
A permutation test of r is equivalent to a permutation test of the OLS slope
P-perm for SMA = NA because the SMA slope cannot be tested
Regression results
Method Intercept Slope Angle (degrees) P-perm (1-tailed)
1 OLS -0.04822033 0.2012451 11.37851 0.01
2 MA -0.18015300 0.2914863 16.25068 0.01
3 SMA -0.64119444 0.6068361 31.25089 NA
Confidence intervals
Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
1 OLS -0.2927963 0.1963556 0.03510125 0.3673889
2 MA -0.5665849 0.1582219 0.06003974 0.5558036
3 SMA -0.9167477 -0.4309431 0.46302535 0.7953130
Eigenvalues: 0.03192833 0.009336979
H statistic used for computing C.I. of MA: 0.04919527 fit3.setosa <- deming(Petal.Width~Petal.Length, setosa)
fit3.setosa
Call:
deming(formula = Petal.Width ~ Petal.Length, data = setosa)
n= 50
Coef se(coef) lower 0.95 upper 0.95
Intercept -0.1801507 0.2016353 -0.57534872 0.2150472
Slope 0.2914848 0.1423112 0.01255992 0.5704096
Scale= 0.0976294 deming_residuals <- resid(fit3.setosa)
model <- fit3.setosa$model
deming_model <- data.frame(model, deming_residuals)
ggplot(subset(setosa), aes(x=Petal.Length, y=deming_residuals))+
geom_point()
The residues look randomly scattered and can be assumed to be normal. There is no calculated R-squared value for the model. You can plot the model against the OLS and SMA models. Here the Deming model is shown in Turquoise and it is closer to the OLS model than the SMA model.
ggplot(subset(setosa), aes(x=Petal.Length, y=Petal.Width))+
geom_point()+
geom_smooth(method=lm, colour="#4b0082")+
stat_ma_line(method="SMA",colour="darkorange") +
geom_abline(slope=0.291, intercept=-0.18, colour="darkturquoise", linewidth=1)