# R/helpers.R In multxpert: Common Multiple Testing Procedures and Gatekeeping Procedures

#### Defines functions bonfsimesincsimespvaltrunc

```# Bonf function is a secondary function which computes the weighted Bonferroni
# p-value for an intersection hypothesis
bonf<-function(p,w)
# P, Vector of raw p-values
# W, Vector of hypothesis weights
{
# Number of null hypotheses included in the intersection
k<-length(w[w!=0])
if (k>0) bonf<-min(p[w!=0]/w[w!=0])
if (k==0) bonf<-1
return(bonf)
}
# End of bonf

# Simes function is a secondary function which computes a truncated version of the
# weighted Simes p-value for an intersection hypothesis
simes<-function(p,w,gamma)
# P, Vector of raw p-values
# W, Vector of hypothesis weights
# GAMMA, Truncation parameter
{
# Number of null hypotheses included in the intersection
k<-length(w[w!=0])
if (k>1)
{
# P-values
temp[1,]<-p[w!=0]
# Weights
temp[2,]<-w[w!=0]
# Normalized weights
temp[3,]<-w[w!=0]/sum(w[w!=0])
# Sort by p-values
sorted<-temp[,order(temp[1,])]
numer<-sorted[1,]
denom<-gamma*cumsum(sorted[3,])+(1-gamma)*sorted[2,]
simes<-min(numer/denom)
}
if (k==1)
{
numer<-p[w!=0]
denom<-gamma+(1-gamma)*w[w!=0]
simes<-numer/denom
}
if (k==0) simes<-1
return(simes)
}
# End of simes

# IncSimes function is a secondary function which computes a truncated version of the weighted
# incomplete Simes p-value for an intersection hypothesis
incsimes<-function(p,w,gamma)
# P, Vector of raw p-values
# W, Vector of hypothesis weights
# GAMMA, Truncation parameter
{
# Number of null hypotheses included in the intersection
k<-length(w[w!=0])
if (k>1)
{
temp<-matrix(0,3,k)
# P-values
temp[1,]<-p[w!=0]
# Weights
temp[2,]<-w[w!=0]
# Normalized weights
temp[3,]<-w[w!=0]/sum(w[w!=0])
# Sort by p-values
sorted<-temp[,order(temp[1,])]
modw<-w[w!=0]
modw[1]<-0
modw[2:k]<-sorted[3,1:k-1]
numer<-sorted[1,]
denom<-gamma*sorted[3,]/(1-cumsum(modw))+(1-gamma)*sorted[2,]
incsimes<-min(numer/denom)
}
if (k==1)
{
numer<-p[w!=0]
denom<-gamma+(1-gamma)*w[w!=0]
incsimes<-numer/denom
}
if (k==0) incsimes<-1
return(incsimes)
}
# End of incsimes

# PValTrunc function is a secondary function which computes adjusted p-values for
# truncated versions of the Holm, Hommel, Hochberg, fixed-sequence and fallback procedures
# in general hypothesis testing problems (equally or unequally weighted null hypotheses)
pvaltrunc<-function(rawp,weight,proc,gamma)
# RAWP, Vector of raw p-values
# WEIGHT, Vector of hypothesis weights
# PROC, Procedure name
# GAMMA, Truncation parameter
{
# Number of null hypotheses
nhyps<-length(rawp)

# Number of weights
nweis<-length(weight)

# Number of intersection hypotheses
nints<-2**nhyps-1

# Decision matrix
h<-matrix(0,nints,nhyps)
for(i in 1:nhyps)
{
for(j in 0:(nints-1))
{
k<-floor(j/2**(nhyps-i))
if (k/2==floor(k/2)) h[j+1,i]<-1
}
}

# Temporary vector of weights
tempw<-matrix(0,1,nhyps)

# Holm procedure
if (proc=="Holm")
{
# Local p-values for Holm procedure
holmp<-matrix(0,nints,nhyps)

for(i in 1:nints)
{
tempw<-(weight*h[i,]/sum(weight*h[i,]))*gamma+(1-gamma)*weight*h[i,]
holmp[i,]<-h[i,]*bonf(rawp,tempw)
}

holm<-rep(0,nhyps)
for(i in 1:nhyps) holm[i]<-pmin(1, max(holmp[,i]))

# Delete temporary objects
rm(holmp)

res<-holm
}

# Hommel procedure
else if (proc=="Hommel")
{
# Local p-values for Hommel procedure
hommp<-matrix(0,nints,nhyps)

for(i in 1:nints) hommp[i,]<-h[i,]*simes(rawp,weight*h[i,],gamma)

hommel<-rep(0,nhyps)
for(i in 1:nhyps) hommel[i]<-max(hommp[,i])

# Delete temporary objects
rm(hommp)

res<-hommel
}

# Hochberg procedure
else if (proc=="Hochberg")
{
# Local p-values for Hochberg procedure
hochp<-matrix(0,nints,nhyps)

for(i in 1:nints) hochp[i,]<-h[i,]*incsimes(rawp,weight*h[i,],gamma)

hochberg<-rep(0,nhyps)
for(i in 1:nhyps) hochberg[i]<-max(hochp[,i])

# Delete temporary objects
rm(hochp)

res<-hochberg
}

# Fallback procedure
else if (proc=="Fallback")
{
# Local p-values for fallback procedure
fallp<-matrix(0,nints,nhyps)
fallwin<-matrix(0,nhyps,nhyps)

# Compute window variables
for(i in 1:nhyps)
{
for(j in 1:nhyps)
{
if (i>=j) fallwin[i,j]<-1
}
}

for(i in 1:nints)
{
tempw[]<-0
for(j in 1:nhyps)
{
if (h[i,j]==1) tempw[j]<-(sum(weight*fallwin[j,])-sum(tempw))
if (h[i,j]==0) tempw[j]<-0
}
tempw<-tempw*gamma+(1-gamma)*weight*h[i,]
fallp[i,]<-h[i,]*bonf(rawp,tempw)
}

fallback<-rep(0,nhyps)
for(i in 1:nhyps) fallback[i]<-pmin(1,max(fallp[,i]))

# Delete temporary objects
rm(fallp)
rm(fallwin)

res<-fallback
}

# Fixed-sequence procedure
else if (proc=="Fixed-sequence")
{
# Compute adjusted p-values using a recursive algorithm
fixedseq<-rep(0,nhyps)
fixedseq[1]<-rawp[1]
if (nhyps>1)
{
for(i in 2:nhyps) fixedseq[i]<-max(fixedseq[i-1],rawp[i])
}

res<-fixedseq
}

rm(tempw)
return(res)
}
# End of pvaltrunc
```

## Try the multxpert package in your browser

Any scripts or data that you put into this service are public.

multxpert documentation built on Oct. 2, 2019, 5:02 p.m.