Monday, October 25, 2010

R: SMA (standard major axis) Regression

# !!! One should first test the significance of the correlation coefficient (r) to determine if the hypothesis of a relationship is supported. No SMA regression equation should be computed when this condition is not met. This remains a less-than-ideal solution since SMA slope estimates cannot be tested for significance. Confidence intervals should also be used with caution: simulations have shown that, as the slope departs from +/-1, the SMA slope estimate is increasingly biased and the confidence interval includes the true value less and less often. Even when the slope is near +/-1, the confidence interval is too narrow if n is very small or if the correlation is weak !!!

# Data input: generic

# To test the significance of the correlation coefficient (r):
cor.test(infile$var1,infile$var2,method="spearman")

# If the result of the above is significant(ie: your variables are significantly correlated), then the following function may be applied.
# ref for regression method: http://www.palass.org/modules.php?name=palaeo_math&page=7

sma.reg <- function(X,Y) # where X is the independent variable and Y is the dependent variable.
{
# calculate the regression slope:
std.devX <- sd(X)
std.devY <- sd(Y)
reg.slope <- std.devX/std.devY

# calculate the y-intercept
av.X <- mean(X)
av.Y <- mean(Y)
y.intercept <- av.Y-(reg.slope*av.X)

# calculate the sign of the slope:
Xmc <- X-av.X
Ymc <- Y-av.Y
XmcYmc <- Xmc*Ymc
slope.dir <- sum(XmcYmc)
if(slope.dir<0)
{
min.Ycoords <- -1*reg.slope*(min(X))+y.intercept
max.Ycoords <- -1*reg.slope*(max(X))+y.intercept
Xcoords <- c(min(X),max(X))
Ycoords <- c(min.Ycoords,max.Ycoords)
plot(X,Y)
lines(Xcoords,Ycoords)
paste("Y","=","-",reg.slope,"X","+",y.intercept)
} else
{
min.Ycoords <- reg.slope*(min(X))+y.intercept
max.Ycoords <- reg.slope*(max(X))+y.intercept
Xcoords <- c(min(X),max(X))
Ycoords <- c(min.Ycoords,max.Ycoords)
plot(X,Y)
lines(Xcoords,Ycoords)
paste("Y","=",reg.slope,"X","+",y.intercept)
}
}
# /* END OF SMA REGRESSION */

No comments:

Post a Comment

Post a Comment