Nothing
##########################################################################################
#Under normality assumption
#Omnibus Chi-square test compare >=3 markers to examine if all have statistically equivalent summary measures
#When exactly 2 markers are compared, z-test on equivalence in summary measures between two markers
#
#########################################################################################
DiagTest3Grp.Test <- function(dat,paired=FALSE,type=c("VUS","Youden"),p=0,q=0,mu=0,conf.level=0.95,alternative=c("two.sided","less","greater"))
{
###Inputs:
#(1) dat: For unpaired data (see paired argument), dat is a list of K components, each corresponding to a marker. Each component of the list should be a data frame of n observation by 2 columns: group must be at 1st column with levels of D-, D0 and D+ corresponding to increasing disease severity and numeric marker measurements at 2nd column ; for paired data, dat should be a data frame with observations at row and group/markers at column (group must be the 1st column with levels D-, D0. D+, followed by markers at 2nd to the last column
#(2)paired: whether all markers are measured on the same set (paired=TRUE) or not (paired=FALSE)
#(3)type: which summary measure (VUS or Youden) will be used to indicate the discriminative ability of the three-group diagnostic tests. Allows unique argument matching,i.e., can use "V" , "Y" as abbreviation
#(4)p, q: if testing on "VUS" and partial VUS satisfying specificity >=p and sensitivity >=q is desired, default=0 for full VUS.
#(5):conf.level: default 95% confidence interval on point estimate (abs diff between two markers's summary measures)
##(6)alternative:fixed as "two.sided" for testing >=3 markers.For testing two markers, can be "two.sided/less/greater" (marker1 at first column-marker2 at second column)
###outputs:
###(1)statistic: for independent samples, statistic is the calculated Chi-square-statistic with df=K-1;for paired samples, statistic is the calculated z-statistic following N(0,1)
###(2)p.value: calculated p-value
###(3)estimate: the markers' summary measures
###(4)conf.int: the 95% CI on the absoluate difference between the two markers' summary measures, only returned for two marker (independent samples or paired samples) , for >=3 markers, NA
###(5)method: "VUS" or "Youden", which summary measure was tested
###(6)
BetweenMarker.cov <- function(marker1, marker2,n.minus,n0,n.plus,rho.minus,rho0,rho.plus,type="VUS")
{
###Calculate the covariance between the two markers' summary measures: cov(theta_1 and theta_2)
##input: (1)marker1, marker2: each a data frame (for marker 1, 2 respectively) containing: a, b, c, d,(mostly a column of res0$patialDeriv outputted above) and V.a, V.b, V.c, V.d
## (2)n.minus, n0,n.plus: the number of obs in each group that both markers are available
#calculate the i<=4, j>4 between marker terms for eq (8) on P686 of Xiong et al on Biometrical,
###output: covariance
mu.s.cov <- function(n,rho,s1,s2)
{
##the covariance between the the two normal mean parameters in the same Diagnosis group for marker 1&2, e.g. in the D- group, cov(mu_{-1},mu_{-2})=rho_{-12}*s_{-1}*s{-2}/n_{-}
mu.cov <- rho*s1*s2/n
##the covariance between the mu parameter in the same Diagnosis group for marker 1&2,e.g. in the D- group, cov(s_{-1},s_{-2})=s_{-1}*s{-2}*rho_{-12}^2*(1-rho_{-12}^2)^3 / [2*n*(3*rho_{-12}^4-rho_{-12}^6-3*rho_{-12}^2+1)]
s.cov <- s1*s2*rho^2*(1-rho^2)^3 /(2*n*(3*rho^4-rho^6-3*rho^2+1))
return(c(mu.cov=mu.cov,s.cov=s.cov))
}
if(type=="VUS")
{
##the same param but two markers (index 1, 2 for marker 1, 2)
a1.a2.cov <- marker1$a*marker2$a*(rho0^2/(2*n0)+rho.minus^2/(2*n.minus))
b1.b2.cov <- rho.minus^2*marker1$b*marker2$b/(2*n.minus)+rho0*marker1$a*marker2$a/n0+rho.minus/n.minus##
c1.c2.cov <- marker1$c*marker2$c*(rho0^2/(2*n0)+rho.plus^2/(2*n.plus))
d1.d2.cov <- rho.plus^2*marker1$d*marker2$d/(2*n.plus)+rho0*marker1$c*marker2$c/n0+rho.plus/n.plus##
##diff para
a1.b2.cov <- rho.minus^2*marker1$a*marker2$b/(2*n.minus)
a1.c2.cov <- rho0^2*marker1$a*marker2$c/(2*n0)
b1.d2.cov <- rho0*marker1$a*marker2$c/n0
c1.d2.cov <- rho.plus^2*marker1$c*marker2$d/(2*n.plus)
###
b1.a2.cov <- rho.minus^2*marker1$b*marker2$a/(2*n.minus)
c1.a2.cov <- rho0^2*marker1$c*marker2$a/(2*n0)
d1.b2.cov <- rho0*marker1$c*marker2$a/n0
d1.c2.cov <- rho.plus^2*marker1$d*marker2$c/(2*n.plus)
cov0 <- marker1$V.a*marker2$V.a*a1.a2.cov+marker1$V.b*marker2$V.b*b1.b2.cov+marker1$V.c*marker2$V.c*c1.c2.cov+marker1$V.d*marker2$V.d*d1.d2.cov+marker1$V.a*marker2$V.b*a1.b2.cov+marker1$V.a*marker2$V.c*a1.c2.cov+marker1$V.b*marker2$V.d*b1.d2.cov+marker1$V.c*marker2$V.d*c1.d2.cov +marker1$V.b*marker2$V.a*b1.a2.cov+marker1$V.c*marker2$V.a*c1.a2.cov+marker1$V.d*marker2$V.b*d1.b2.cov+marker1$V.d*marker2$V.c*d1.c2.cov
}
else
{
res.minus <- mu.s.cov(n.minus,rho.minus,marker1$s.minus,marker2$s.minus)
res0 <- mu.s.cov(n0,rho0,marker1$s0,marker2$s0)
res.plus <- mu.s.cov(n.plus,rho.plus,marker1$s.plus,marker2$s.plus)
mu.minus.cov <- res.minus["mu.cov"]
s.minus.cov <- res.minus["s.cov"]
mu0.cov <- res0["mu.cov"]
s0.cov <- res0["s.cov"]
mu.plus.cov <- res.plus["mu.cov"]
s.plus.cov <- res.plus["s.cov"]
cov0 <- marker1$Y.mu.minus*marker2$Y.mu.minus*mu.minus.cov+marker1$Y.mu0*marker2$Y.mu0*mu0.cov+marker1$Y.mu.plus*marker2$Y.mu.plus*mu.plus.cov+marker1$Y.s.minus*marker2$Y.s.minus*s.minus.cov+marker1$Y.s0*marker2$Y.s0*s0.cov+marker1$Y.s.plus*marker2$Y.s.plus*s.plus.cov
}
return(cov0)
}
##################Start ###################
type <- match.arg(type)##allows incomplete input in the argument, e.g., "V" for "VUS" and "Y" for "Youden"
if(!type%in%c("VUS","Youden")) stop("type must be VUS or Youden!!")
dat.str <- class(dat)
if(dat.str=="list" & paired) stop("Paired data (where markers measured on the same setof subjects) should be organized as a data frame with samples at row and group& marker measurements at column, not a list!")
if(dat.str=="data.frame" & !paired) stop("unpaired data (where markers measured on different sets subjects) should be organized as a list (with each component for a marker), not a data frame!")
####check data structure
if(dat.str=="list") K <- length(dat)
if(dat.str=="data.frame") K <- ncol(dat)-1##the first column is group membership
if(K<2) stop("must have >=2 markers!")
alternative <- match.arg(alternative)
data.name <- paste("Test ", deparse(substitute(dat))," on ", type,sep="")
###Construct the contrast matrix of K*(K-1) at row r, the r-th position is 1 and the (r+1) position is -1
A <- matrix(0,K-1,K)
diag(A) <- 1
for(ii in 1:(K-1)) A[ii,ii+1] <- -1
###independent samples
if(!paired)
{
###Obtain point estimate on the summary measure and its associated variance
markerID <- names(dat)
res0 <- lapply(dat,function(marker.dat)
{
marker.list <- split(marker.dat[,2],marker.dat[,1])##make sure group at 1st column and marker measurements at 2nd column
xx <- na.exclude(marker.list$'D-')
yy <- na.exclude(marker.list$D0)
zz <- na.exclude(marker.list$'D+')
if(type=="VUS")
{
vus.res <- VUS(x=xx,y=yy,z=zz,method="Normal",p=p,q=q)
data.frame(estimate=vus.res$estimate,variance=vus.res$variance,vus.res$partialDeriv)
}
else
{
youden.res <- Youden3Grp(x=xx,y=yy,z=zz,method="Normal")
data.frame(estimate=youden.res$estimate,variance=youden.res$variance,youden.res$partialDeriv)
}
})
#browser()
res0 <- do.call("rbind",res0)
theta <- res0$estimate
variance <- res0$variance
###construct the covariance matrix on theta estimates with variance on estimate at the diagonal and all other entries are 0
Sigma <- diag(variance)
}
###paired samples, need to calculate the covariance
if(paired)
{
markerID <- names(dat)[-1]
###Obtain point estimate on the summary measure and its associated variance
res0 <- apply(dat[,-1],2,function(marker.dat)
{
marker.list <- split(marker.dat,dat[,1])##make sure group at 1st column and marker measurements at 2nd column
xx <- na.exclude(marker.list$'D-')
yy <- na.exclude(marker.list$D0)
zz <- na.exclude(marker.list$'D+')
if(type=="VUS")
{
vus.res <- VUS(x=xx,y=yy,z=zz,method="Normal",p=p,q=q)
#partialDeriv <- vus.res$partialDeriv
data.frame(estimate=vus.res$estimate,variance=vus.res$variance,vus.res$partialDeriv)
}
else
{
youden.res <- Youden3Grp(x=xx,y=yy,z=zz,method="Normal")
data.frame(estimate=youden.res$estimate,variance=youden.res$variance,youden.res$partialDeriv)
}
})
#browser()##make sure the output is a data frame with the names
res0 <- do.call("rbind",res0)
theta <- res0$estimate
variance <- res0$variance
partialDeriv <- subset(res0,select=-c(estimate,variance))##a data frame of K row (each corresponding to a marker) and columns are (1) partial derivative w.r.t a, b, c, d: V.a, V.b, V.c, V.d (2) variance, covariance: a.var, b.var,c.var,d.var,a.b.cov,a.c.cov,b.d.cov,c.d.cov
###obtain within each diagnosis group correlation among markers
cor.minus <- cor(dat[dat[,1]=="D-",-1],use="pairwise.complete.obs")##a matrix of K*K, cor within D- group
cor0 <- cor(dat[dat[,1]=="D0",-1],use="pairwise.complete.obs")
cor.plus <- cor(dat[dat[,1]=="D+",-1],use="pairwise.complete.obs")
###construct the covariance matrix by placing the variance on theta estimates at the diagnonal term
Sigma <- diag(variance)
####put at the entry (i,j) of Sigma the covariance of marker i and marker j
if(K==2)
{
###get sample size with complete obs on each pair of marker and have group label
temp.dat <- na.exclude(dat[,c(1,2,3)])
n.minus <- sum(temp.dat[,1]=="D-")
n0 <- sum(temp.dat[,1]=="D0")
n.plus <- sum(temp.dat[,1]=="D+")
#browser()
current.cov <- BetweenMarker.cov(marker1=partialDeriv[1,], marker2=partialDeriv[2,],n.minus=n.minus,n0=n0,n.plus=n.plus,rho.minus=cor.minus[1,2],rho0=cor0[1,2],rho.plus=cor.plus[1,2],type=type)
Sigma[1,2] <- Sigma[2,1] <- current.cov
}else{
for (i1 in 1:(K-1))
{
for(i2 in (i1+1):K)
{
temp.dat <- na.exclude(dat[,c(1,i1,i2)])
n.minus <- sum(temp.dat[,1]=="D-")
n0 <- sum(temp.dat[,1]=="D0")
n.plus <- sum(temp.dat[,1]=="D+")
current.cov <- BetweenMarker.cov(marker1=partialDeriv[i1,], marker2=partialDeriv[i2,],n.minus=n.minus,n0=n0,n.plus=n.plus,rho.minus=cor.minus[i1,i2],rho0=cor0[i1,i2],rho.plus=cor.plus[i1,i2],type=type)
Sigma[i1,i2] <- Sigma[i2,i1] <- current.cov
}
}
}
}##end of paired
#print(Sigma)
###calculate A*theta
A.theta <- A%*%matrix(theta,ncol=1)
###Calcualte the Chi-square stat and associated P-value
A.Sigma.A <- A%*%Sigma%*%t(A)
inv.mat <- solve(A.Sigma.A)
#print(inv.mat)
chisq.stat <- t(A.theta)%*%inv.mat%*%A.theta
names(chisq.stat) <- "Chi-square"
df <- K-1
names(df) <- "df"
chisq.p.value <- pchisq(chisq.stat,df=K-1,lower.tail=FALSE)
names(mu) <- paste("diff in ",type,sep="")
####For two markers, do z-test
estimate <- theta
names(estimate) <- paste(paste(type," of ",sep=""),markerID,sep="")
if(K==2)
{
diff <- theta[1]-theta[2]
diff.se <- sqrt(as.numeric(A.Sigma.A))
###alternative=="two.sided"
z.stat <- (diff-mu)/diff.se
names(z.stat) <- "Z-stat"
if(alternative=="two.sided")
{
z.p.value <- 2*pnorm(-abs(z.stat),lower.tail=TRUE)
#ci <- CI.normal(diff,abs.diff.se,alpha=1-conf.level)
alpha <- 1-conf.level
cint <- qnorm(1-alpha/2)*diff.se
cint <- diff+c(-cint,cint)
}
else if (alternative=="less")
{
z.p.value <- pnorm(z.stat,lower.tail=TRUE)
cint <- c(-Inf,diff+qnorm(conf.level)*diff.se)
}else if (alternative=="greater")
{
z.p.value <- pnorm(z.stat,lower.tail=FALSE)
cint <- c(diff-qnorm(conf.level)*diff.se,Inf)
}
attr(cint,"conf.level") <- conf.level
method <- "normal-test"
mean.sd <- c(mean=0,sd=1)
rval <- list(statistic=z.stat,p.value=z.p.value,parameter=mean.sd,estimate=estimate,conf.int=cint,method=method,type=type,null.value=mu,alternative=alternative,data.name=data.name,Sigma=Sigma)
}else{
method <- "Chi-Sqaure test"
rval <- list(statistic=chisq.stat,p.value=chisq.p.value,parameter=df,estimate=estimate,conf.int=NA,method=method,type=type,null.value=mu,alternative=alternative,data.name=data.name,Sigma=Sigma)
}
class(rval) <- "htest"
return(rval)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.