#' MULTIVARIATE PERMUTATION TESTING OF CORRELATIONS
#'
#' "This R program performs multivariate permutation tests of the correlation between a single criterion variable and multiple predictor variables. Data must be imported by user with subjects/cases as rows and variables as columns. Missing data must be specified."
#' This script was Written by Jennifer Urbano Blackford, Ph.D. of The Kennedy Center, Vanderbilt University, Nashville, TN.
#' The original can be found at http://users.cla.umn.edu/~nwaller/downloads/mpt/mptcorr.r
#'
#' @export
#----------------------------------------------------------------------------
#MULTIVARIATE PERMUTATION TESTING OF CORRELATIONS
#Description: This R program performs multivariate permutation tests of
# the correlation between a single criterion variable and
# multiple predictor variables.
# Data must be imported by user with subjects/cases as rows
# and variables as columns. Missing data must be specified.
#
#Written by: Jennifer Urbano Blackford, Ph.D.
#
#Owned by: The Kennedy Center, Vanderbilt University, Nashville, TN
#--------------------------------------------------------------------------
#define function & specify arguments
#defaults are specified with an '='
MPT.Corr<-function(observed.data,criterion,predictor.start,predictor.end,
output.file="Results.txt",test=2,tail=2,alpha=.05,output=0,permutations=10000){
#--------------------------------------------------------------------------
#CREATE NECESSARY DATASETS AND VARIABLES
#set user specifications
#set text to print to screen or print to specified file
if (output.file != '') (sink(output.file))
#create datasets using the user specifications above
predict.data<-as.matrix(observed.data[, predictor.start:predictor.end,drop=FALSE])
if (output==1) (print (predict.data))
crit.data<-as.matrix(observed.data[,criterion,drop=FALSE])
if (output==1) (print (crit.data))
#Assign variable numbers from predictor and criterion datasets
#set number of subjects as number of rows in Data
nsubj<-nrow(crit.data)
#set number of predictors to number of columns in predict.data
npredict<-ncol(predict.data)
#set number of criterion variables to number of columns in crit.data
ncrit<-ncol(crit.data)
#save the variable names for later use
colnames <- dimnames(predict.data)[[2]]
#set tail information for the cor.test function
if (test==2) (alt<-"t")
if (test==1 && tail==1) (alt<-"g")
if (test==1 && tail==-1) (alt<-"l")
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#CREATE OBSERVED AND SIMULATED CORRELATION DATASETS
#calculate correlation for actual data--crit.data and predict.data
#save correlations in one matrix and p values in another
observed.r<-matrix(0,1,npredict,dimnames=list(c("r"),c(colnames)))
observed.p<-matrix(0,1,npredict,dimnames=list(c("p"),c(colnames)))
#create a matrix to save the simulated r values
simulated.r<-matrix(0, permutations,npredict)
#create a vector OriginalCorr to save correlation results
OriginalCorr<-apply(predict.data,2,cor.test,y=crit.data,alternative=alt)
for(i in 1:npredict){
observed.r[i]<-as.matrix(OriginalCorr[[i]]$estimate) #get the r value
observed.p[i]<-as.matrix(OriginalCorr[[i]]$p.value)} #get the p value
#print correlation results
writeLines("Multivariate Permutation Test for Correlations")
writeLines(" Yoder, Blackford, Waller, and Kim (2003)")
writeLines("")
writeLines("Observed Correlations & P-values")
writeLines("")
print(cbind(t(round(observed.r,digits=3)),t(round(observed.p,digits=3))))
writeLines("")
#create simulated r by shuffling the original criterion data and computing correlations
for(loop in 1:permutations){
#randomly sampling nsubj times without replacement=shuffling
ShuffleData<-as.matrix(crit.data[sample(1:nsubj),])
#save the post-shuffle corr for each variable and each loop into a matrix called simulated.r
#note: when using the corr function you must specify to use pairwise deletion
SimulatedCorr<-cor(predict.data,ShuffleData,use="pairwise.complete.obs")
for(i in 1:npredict){
simulated.r[loop,i]<-SimulatedCorr[i]}
} #end bracket for(loop in 1:permutations)
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#PERFORM MULTIVARIATE PERMUATION TESTS
#create a matrix to store the MPT results
Results<-matrix(,1,4,dimnames=list(NULL,c("Variable","r","p","MPT exact p")))[-1,]
#create matrices of maximum corrs and observed corrs in each permutation
#loop to test significance for max r (i.e., step down procedure).
while (npredict > 0) {
#select largest of the absolute values--save column position
#using column as placeholder provides a way to track the sign of the r
max.col<-matrix(0,permutations)
for (loop in 1:permutations){
max.col[loop,]<-(1:npredict)[abs(simulated.r[loop,])== max(abs(simulated.r[loop,]))]}
#save the largest value--maintain original sign (using max.col)
max.simulated.r<-matrix(0,permutations)
for (loop in 1:permutations){
max.simulated.r[loop,1]<-simulated.r[loop,max.col[loop,]]}
#print histogram
if (output==1) (hist(max.simulated.r))
#select the largest correlation from the observed correlation for testing
#use test specification information to determine which correlation to use
#set max.col.r to the appropriate column number
#two tailed test: get largest absolute value
if (test==2) (max.col.r <-(order(abs(observed.r[1,]))[npredict]))
#one-tailed positive test: get largest positive value
if (test==1 && tail==1) (max.col.r<-(order(observed.r[1,])[npredict]))
#one-tailed negative test: get largest negative value
if (test==1 && tail==-1) (max.col.r<-(order(observed.r[1,])[1]))
#use the max.col.r variable to get the name, r, and p-value
#associated with the largest correlation to be tested
largest.observed.r<-observed.r[max.col.r]
smallest.observed.p<-observed.p[max.col.r]
variable<-colnames[max.col.r]
#apply breaks to stop the program if necessary
#if testing the positive tail only then stop the program when the observed r is less than 0
if (test==1 && tail==1 && largest.observed.r <0) break
#if testing the negative tail only then stop the program when the observed r is greater than 0
if (test==1 && tail==-1 && largest.observed.r >0) break
#compare largest correlation to the max r distribution
#determine significance by creating an array of 0s and 1s representing
#whether the simulated corr is larger than the observed corr.
p.dist<-rep(0,permutations)
#for two tailed test
#positive observed correlation--check positive tail & negative tail
if (test==2 && (largest.observed.r >0)) (p.dist[max.simulated.r>= largest.observed.r]<-1)
if (test==2 && (largest.observed.r >0)) (p.dist[max.simulated.r<= largest.observed.r*-1]<-1)
#negative observed correlation--check positive tail & negative tail
if (test==2 && (largest.observed.r <0))(p.dist[max.simulated.r<= largest.observed.r]<-1)
if (test==2 && (largest.observed.r <0))(p.dist[max.simulated.r>= largest.observed.r*-1]<-1)
#for one tailed positive test check only the positive tail
if (test==1 && tail==1) (p.dist[max.simulated.r>= largest.observed.r]<-1)
#for one tailed negative test check only the negative tail
if (test==1 && tail==-1) (p.dist[max.simulated.r<= largest.observed.r]<-1)
#get p-value as number of permutations above/below observed r and divided by number of permutations
p.value<-0
if (sum(p.dist) > 0) (p.value<-(sum(p.dist)/permutations))
#stop the program when p.value is greater than the value specified
if (p.value > alpha) break
#if the p-value was significant then delete the column for max r and continue
#delete the column associated with the maximum observed corr
simulated.r<-simulated.r[,-max.col.r,drop=FALSE]
colnames<-colnames[-max.col.r,drop=FALSE]
observed.r<-observed.r[,-max.col.r,drop=FALSE]
observed.p<-observed.p[,-max.col.r,drop=FALSE]
#reset number of variables
npredict<-npredict-1
#create a matrix to store the results
ResultsTemp<-cbind(variable,round(largest.observed.r,digits=3),
round(smallest.observed.p,digits=3),round(p.value, digits=3))
Results<-rbind(Results,ResultsTemp)
rm(ResultsTemp)
#end testing
}
#add row labels to Results
if (nrow(Results)>=1) (dimnames(Results)[[1]]<-c(1:nrow(Results)))
#print final results if there are any
writeLines("Results of MPT for Correlations")
writeLines("")
if (nrow(Results)>=1) print(Results,quote=FALSE,rowlab=c(1:nrow(Results)))
if (nrow(Results)<1) writeLines("No significant Results")
writeLines("")
sink() #close sink
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.