RMA <- function(x, y)
{
CL95intercept <- (intercept <- mean(y)-(slope <- sign(cor(x, y))*sqrt(var(y)/var(x)))*mean(x))+c(1, -1)*qt(0.025, df <- (n1 <- (n <- length(x))-1)-1)*(SEintercept <- sqrt((MSE <- (var(y)-cov(x, y)^2/var(x))*n1/df)*(1/n+mean(x)^2/var(x)/n1)))
CL95slope <- slope+c(1, -1)*qt(0.025, df)*(SEslope <- sqrt(MSE/var(x)/n1))
result1 <- c(slope, SEslope, CL95slope)
result2 <- c(intercept, SEintercept, CL95intercept)
names(result1) <- c("slope", "SE[slope]", "95% LCL", "95% UCL")
names(result2) <- c("intercept", "SE[int.]", "95% LCL", "95% UCL")
list(RMAslope=result1, RMAintercept=result2)
}
y <- c(61, 37, 65, 69, 54, 93, 87, 89, 100, 90, 97)
x <- c(14, 17, 24, 25, 27, 33, 34, 37, 40, 41, 42)
RMA(x,y)
# $RMAslope
# slope SE[slope] 95% LCL 95% UCL
# 2.1193664 0.3324959 1.3672085 2.8715243
# $RMAintercept
# intercept SE[int.] 95% LCL 95% UCL
# 12.19378 10.54975 -11.67141 36.05898