#' VCO2sim Determines the Probability that Two VCO2 Measures are the Same Between Testing Days
#'
#' This function models the reliability of VCO2 data from a particular system as univariate normal
#' based on the work by Crouter 2006
#' Additional credit goes to Stackexchange user Wolfgang for his example code of the overlapping coefficient
#'
#'
#' @import ggplot2 msm cubature
#' @param a The first VCO2 value obtained
#' @param b The second VCO2 value obtained
#' @param system_a The system used, in quotes, to obtain measurement 'a'.
#' Accepts "parvo_2400" and "douglas_bag". Defaults to parvo_2400.
#' @param system_b The system used, in quotes, to obtain measurement 'b'.
#' Accepts "parvo_2400" and "douglas_bag". Defaults to parvo_2400.
#' @param plot False returns the probability they are same distribution. True returns a plot of overlapping distributions.
#' @export
#' @examples
#' put something here
VCO2_sim<- function(a, b, system_a='parvo_2400', system_b='parvo_2400', plot = F){
library(ggplot2)
library(msm)
library(cubature)
mu1 <- a
mu2 <- b
if (system_a == 'parvo_2400'){sd1 <- suppressWarnings(predict(VCO2_parvo_pred, data.frame(VCO2_mean = mu1)))}
if (system_b == 'parvo_2400'){sd2 <- suppressWarnings(predict(VCO2_parvo_pred, data.frame(VCO2_mean = mu2)))}
if (system_a == 'douglas_bag'){sd1 <- suppressWarnings(predict(VCO2_douglas_pred, data.frame(VCO2_mean = mu1)))}
if (system_b == 'douglas_bag'){sd2 <- suppressWarnings(predict(VCO2_douglas_pred, data.frame(VCO2_mean = mu2)))}
#throw warning that algorithm is extrapolating beyond original data
if( any(mu1 > 4.19 | mu2 > 4.19) ) warning('Algorithm extrapolating beyond original data, results may be unreliable/unstable')
#create funtion for determining integral
min.f1f2 <- function(x, mu1, mu2, sd1, sd2, lowerb) {
f1 <- dtnorm(x, mean=mu1, sd=sd1, lower=lowerb)
f2 <- dtnorm(x, mean=mu2, sd=sd2, lower=lowerb)
pmin(f1, f2)
}
lowerb <- 0
if (plot == F ){
#just return integral information
areacurve <- pcubature(min.f1f2, 0, 1000, mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2, lowerb=0, fDim = 4)
prob_diff <- areacurve[[1]][1]*100
#integration occasionally returns tiny negative numbers, turn them to tiny positive numbers
if( prob_diff < 0 ){prob_diff = 0.00000000000000000000000000000001}
return(prob_diff)
#This error is known, but can't determine root cause
if( prob_diff == 0 ) stop('Known error in cubature integration, please re-run function to return non-zero answer.')
}else if (plot ==T){
#Return the simulation data for plotting
#this needs to be fixed so all the data will plot correctly, just do it later
xs <- seq(min(mu1 - 3*sd1, mu2 - 3*sd2), max(mu1 + 3*sd1, mu2 + 3*sd2), .001)
f1 <- dtnorm(xs, mean=mu1, sd=sd1, lower=lowerb)
f2 <- dtnorm(xs, mean=mu2, sd=sd2, lower=lowerb)
ys <- min.f1f2(xs, mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2, lowerb = 0)
overall <- data.frame(f1,f2,xs, ys)
dat_plot <-ggplot(overall, aes(x= xs, y=f1)) +
geom_line() +
geom_line(aes(x=xs, y=f2)) +
geom_ribbon(aes(ymin =0, ymax=ys)) +
labs(
x= expression('VCO'[2]*' (L/min)'),
y= 'Probability Density'
) +
theme(panel.background = element_rect(fill = 'white'))
return(dat_plot)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.