R/ccp.a.r

In jjlynch2/AnthroStats: Anthropology Statistics

#' Derive a p-value for a case comparison in a simple linear regression model
#'  
#' @param X The independent variable from case comparison
#' @param Y The dependent variable from case comparison
#' @param reference_data The reference data for the linear model
#' @param prediction_interval The prediction interval to use for the case comparison
#' @param tails The number of tails to derive a p-value from the t-distribution
#'
#' @keywords ccp.a
#' @export
#' ccp.a()

ccp.a <- function(X = NULL, Y = NULL, reference_data = NULL, prediction_interval = 0.95, tails = 2) {

	if(ncol(reference_data) > 2) {
		print("Only two columns are supported")
		return()
	}
	
	reference_data <- as.data.frame(reference_data)
	nref <- nrow(reference_data)
	colnames(reference_data) <- c("A","B")

	lmodel <- lm(A~B, data = reference_data)


	pm1 <- predict(lmodel, newdata = data.frame(B = X), interval = "prediction", level = prediction_interval)

	tt <- abs(pm1[1] - Y) / ( summary.lm((lmodel))$sigma * sqrt( 1+(1/nref) + ((X - mean(reference_data[,1]))^2) / (nref * sd(reference_data[,1])^2) ) )
	pp <- tails * pt(-abs(tt), df = nref - 2)

	a <- data.frame(tt, pp)
	colnames(a) <- c("t-statistic", "p-value")
	return(list(lmodel, a))
}