While at the Joint Statistical Meeting a few weeks ago I was talking to a friend about various aspects to clinical trials. He indicated that no current R package was able to perfectly reproduce Passing-Bablok (PB) regression so that it exactly matched SAS. He ultimately wrote a couple of functions and kindly shared them with me so that I could share them here.
There is currently an R package (mcr) that will do the Passing-Bablok regression but, as it turns out, the output does not precisely match SAS. Sadly, I don’t currently have a copy of SAS available to use for this purpose so I can’t independently run code and comparisons. Though I would be interested in an equivalent comparison. However, I ran some examples below to compare the R mcr package.
These code snippets that he shared with me will only calculate the slope, intercept and the confidence interval for each coefficient. It’s not CRAN package ready but the short function and output are convenient and easy to use.
One interesting characteristic comparing the two PB regressions (at least in the simulation below) is that the mcr package always has a slightly smaller coefficient. In the simulation below the maximum difference is .00054. The minimum difference is somewhere along the lines of zero (1.59e-12). The details of the differences is something I must look into at a later date.
library(mcr)
# Function to generate some correlated data
genData = function(numobs=100,randval){
R = matrix(cbind(1,.80, .80,1),nrow=2)
U = t(chol(R))
nvars = dim(U)[1]
set.seed(randval)
random.normal = matrix(rnorm(nvars*numobs,100,10), nrow=nvars, ncol=numobs);
X = U %*% random.normal
newX = t(X)
raw = newX
return(raw)
}
PB.reg = function(X,Y,alpha=.05) {
## (1) calculate all pairwise slopes
x < - X
y <- Y
dat <- cbind(x,y)
n <- length(x)
S <- array(NA,dim=rep(n,2))
for(i in 1:(n-1)){
for(j in (i+1):n) {
if(i != j) {
S[i,j] <- (y[i] - y[j])/(x[i] - x[j])
}
}
}
S <- sort(na.exclude(as.vector(S)))
K <- sum(S <= -1) - .5 * sum(S == -1)
N <- length(S)
b <- ifelse(N%%2,S[(N+1)/2+K],mean(S[N/2+K+0:1]))
### CI for b
C.gamma <- qnorm(1-alpha/2) * sqrt(n*(n-1)*(2*n+5)/18)
M1 <- round((N-C.gamma)/2,0)
M2 <- N - M1 + 1
CI.b <- c(LB=S[M1+K],UB=S[M2+K])
a <- median(y - b*x)
### CI for a
CI.a <- c(LB=median(y - CI.b["UB"]*x),
UB=median(y - CI.b["LB"]*x))
##
return(list(a=a, CI.a=CI.a, b=b, CI.b=CI.b))
}
# A quick function to give only b coefficient
PB.fit.red = function(X,Y){
fit.lm = PB.reg(X, Y)
fit.lm$b
}
# Generate some one-time use correlated data
simData = genData(numobs=20,c(12345))
# Updated function
pb.fit = PB.reg(simData[,1], simData[,2])
# Look at just the b coefficient, single example
PB.fit.red(simData[,1], simData[,2])
# mcreg function, single example
pb.fit2 = mcreg(simData[,1], simData[,2], method.reg="PaBa", method.ci="bootstrap")
slot(pb.fit2, "glob.coef")[2]
&amp;nbsp;
# Run a small simulation
nsims = 250
sim.data = matrix(NA, nrow=nsims, ncol=2)
# for() loop used for demonstration purposes.
# Not efficient code applying a function over the vectors is more efficient
for(i in 1:nsims){
simData = genData( numobs=20, i )
# Look at just the b coefficient, multiple simulations
sim.data[i,1] = PB.fit.red(simData[,1], simData[,2])
# mcreg function, multiple simulations
pb.fit2 = mcreg(simData[,1], simData[,2], method.reg="PaBa", method.ci="bootstrap")
sim.data[i,2] = slot(pb.fit2, "glob.coef")[2]
}
delta = sim.data[,1] - sim.data[,2]
mean(delta)
max(delta)
min(delta)
hist(delta, nclass=50, main="Distribution of Differences in PB Regression")
plot(sim.data, main="Difference Between Two Types of Passing-Bablok Non-Parametric Regression",
xlab="Improved PB Function", ylab="mcr Package")
[/sourcecode]
) and concluding the alternative (
). The 5% that is often used is known as a Type I error. It is often identified by the Greek letter alpha (
). This 5% is the probability of making a Type I error and is often called significance level. This means that the probability of an error and rejecting the null hypothesis when the null hypothesis is in fact true is 5%.
Follows a Chi Square Distribution
-Variance From two Populations should Be Equal (unless using Welch's test).
-Two Populations Should Be Independent
Two independent populations assumption gets a little sticky unless you truly are looking at two different populations from your data (e.g. Male/Female or Hispanic/Non-Hispanic). In other words just because you're comparing two different questions from your questionnaire doesn't necessarily mean you have two independent populations.
[sourcecode language="css"]
> t.test(raw.1.2$value ~ raw.1.2$group, var.equal=TRUE)
Two Sample t-test
data: raw.1.2$value by raw.1.2$group
t = 2.5622, df = 198, p-value = 0.01114
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.1151745 0.8848255
sample estimates:
mean in group 1 mean in group 2
3.4 2.9
