Wednesday, 3 May 2023

R function slopeRatio

 The slopeRatio function calculates the relative potency of two samples based on their linear regression slopes.

It also returns the confidence interval for the relative potency

The function slopeRatio takes four arguments: x, y0, y1, and alpha

  • x is the independent variable (e.g., dose, time, ...)
  • y0 is the standard sample
  • y1 is the test sample
  • alpha is the significance level for the confidence interval (default is 0.05)


slopeRatio uses the R function lm to fit linear models for y0 and y1 as functions of x. It then extracts the slopes and standard errors of the models and computes the ratio of the slopes (R) and its standard error using a formula from Finney (1978).

It also checks if the denominator slope (B0) is significantly different from zero using a g-value.

If g > 1, then the use of the ratio is not sensible for the data.

It then calculates the lower and upper bounds of the confidence interval for R. 

It returns a list with four elements: Slope, St.Err, RelativePotency, and CI, where:

Slope is a vector of two values, which are the slopes of the regression lines for y0 and y1 against x.

St.Err is a vector of two values, which are the standard errors of the slopes.

RelativePotency is a scalar value, which is the ratio of the slopes of y1 and y0.

CI is a vector of two values, which are the lower and upper bounds of the confidence interval for the relative potency. 


slopeRatio <- function(x,y0,y1,alpha=0.05) {

#y0: Standard sample

#y1: Test sample

#x: independent variable (e.g., dose, time, ...)

N <- length(x)

nY <- 2

DF <- data.frame(cbind(x,y0,y1))

B <- sapply(2:(nY+1), function(i) summary(mod <- lm(DF[,i] ~ x,data=DF))$coef[2]) #slopes of regression lines

B0 <- B[1]

B1 <- B[2]

SE <- sapply(2:(nY+1), function(i) summary(mod <- lm(DF[,i] ~ x,data=DF))$coef[[4]])#St.Err of slopes

SE0 <- SE[1]

SE1 <- SE[2]

R <- B[2]/B[1] #ratio of the slopes

mod0 <- lm(y0 ~ x)

mod1 <- lm(y1 ~ x)

Cov <- abs(summary(lm(y1 ~ y0))$cov.unscaled[1,2])

T <- qt(1-alpha/2,N-1)

g <- (T^2)*(SE0^2)/(B0^2)#if g>1 the denominator, B0, is not significantly different from 0 and the use of the ratio is not sensible for those data

CI.low <- (R-(g*Cov)/(SE0^2)-(T/B0)*sqrt((1-g)*(SE1^2))+(R^2)*(SE0^2)-2*R*Cov +g*(Cov^2)/(SE0^2))/(1-g)

CI.up <- (R-(g*Cov)/(SE0^2)+(T/B0)*sqrt((1-g)*(SE1^2))+(R^2)*(SE0^2)-2*R*Cov +g*(Cov^2)/(SE0^2))/(1-g)

return(list(Slope=B,St.Err=SE,RelativePotency=R,CI=c(CI.low,CI.up)))

}


#EXAMPLE

x<- seq(-6,0,-1)

yS<-c(7.92600, 6.72275, 6.07275, 5.33000, 4.35925, 3.36875, 2.30125)

yT<- c(6.913, 6.018, 5.505, 4.966, 3.440, 3.154, 1.802)

result<-slopeRatio(x,yS,yT)

print(result$Slope)

[1] -0.9034196 -0.8259286

print(result$St.Err)

[1] 0.02979541 0.05703906

print(result$RelativePotency)

[1] 0.9142247

print(result$CI)

[1] -0.6217680 -0.9317596



Since both the bounds of the confidence interval of the relative potency (absolute value) are less than 1, it means that the standard sample is more potent than the test preparation.

References

Grimm, H. (1979), FINNEY, D. J.: Statistical Method in Biological Assay. 3. ed. Charles Griffin & Co., London and High Wycombe 1978. VII, 508 S., £ 19 net. Biom. J., 21: 689-690. https://doi.org/10.1002/bimj.4710210714



No comments:

Post a Comment

Understanding Anaerobic Threshold (VT2) and VO2 Max in Endurance Training

  Introduction: The Science Behind Ventilatory Thresholds Every endurance athlete, whether a long-distance runner, cyclist, or swimmer, st...