Nothing
rand.test <-
function(set1, set2, sims=1000, crit=.95, graph=TRUE, seed=2) {
set1 <- data.frame(set1)
set2 <- data.frame(set2)
# This part gets the data ready to be used in the randomization test
samp.distr=c() #Create a sampling distribution vector for the average absolute r
samp.distsig=c() #Create a sampling distribution vector for the number statistically significant
complete = complete.cases(cbind(set1,set2)) # Combine the data sets and keep only complete cases
set1.set = subset(set1, subset=complete) #Store the "complete" data sets
set2.set = subset(set2, subset=complete)
n = nrow(set1.set) #Find the sample size
critT = qt(.025, n-2, lower.tail=FALSE) # Find the critical t (assumes alpha = .05) for each test
critr = sqrt( critT^2 / (critT^2 + n - 2) ) # Find the critical r value
AbsRObs = mean(abs(cor(set1.set, set2.set))) #Find the Avg. Absolute R Obs
SigObs = sum(abs(cor(set1.set, set2.set)) >= critr) #Find the number significant observed
# This part starts the randomization
if(seed!=F) {set.seed(seed)}
for (i in 1:sims) {
rand.order = sample(n, n, replace=FALSE) #Generate a sample of random orders
cor.mat = cor(set1.set[rand.order,],set2.set) #Get the simulated correlation matrix
samp.distr[i] = mean(abs(cor.mat)) #Store the absolute average simulated r's in samp.distr
samp.distsig[i] = sum(abs(cor.mat) >= critr) #Store the number significant in samp.distsig
}
# This part computes the statistical properties of the two sampling distributions
SimMeanR = mean(samp.distr) #Compute the mean of the sampling distribution
SimSDr = sd(samp.distr) #And the SD
Crit95r = quantile(samp.distr,crit) #And the critical value (default 95th percentile)
pr = sum(samp.distr >= AbsRObs) / sims #Find the probability of the observed value
pr.me <- sqrt(pr * (1-pr) / sims) * qnorm(.9995)
SimMeanSig = mean(samp.distsig) #Compute the mean
SimSDsig = sd(samp.distsig) # SD
Crit95Sig = quantile(samp.distsig,crit) # Critical value (default 95th percentile)
pSig = sum(samp.distsig >= SigObs) / sims # Compute a probability value
pSig.me <- sqrt(pSig * (1-pSig) / sims) * qnorm(.9995)
if(pr + pr.me > 1.00 | pr - pr.me < .00) {warning("Confidence intervals for p-values may be inaccurate. Try a larger number of sims.")}
if(pSig + pSig.me > 1.00 | pSig - pSig.me < .00) {warning("Confidence intervals for p-values may be inaccurate. Try a larger number of sims.")}
if(pr == .00 | pSig == .00) {warning("When p=.00, confidence interval for p not valid.")}
#Clean up and print the results
out.AbsR = round(rbind(n, AbsRObs, SimMeanR, SimSDr, pr, pr + pr.me, pr - pr.me, Crit95r),4)
colnames(out.AbsR) = c("Average Absolute r")
rownames(out.AbsR) = c("N", "Observed", "Exp. By Chance", "Standard Error", "p", "99.9% Upperbound p", "99.9% Lowerbound p", "95th %")
out.Sig = round(rbind(n, SigObs, SimMeanSig, SimSDsig, pSig, pSig + pSig.me, pSig - pSig.me, Crit95Sig),4)
colnames(out.Sig) = c("Number Significant")
rownames(out.Sig) = c("N", "Observed", "Exp. By Chance", "Standard Error", "p", "99.9% Upperbound p", "99.9% Lowerbound p", "95th %")
results <- list("AbsR"=out.AbsR, "Sig"=out.Sig)
# This part creates histogram graphics of the sampling distributions
if (graph == TRUE) {
old.par = par(mfrow=c(2,1)) # Sets the PAR command two produce two vertical histograms
hist(samp.distr, freq=TRUE, col="cyan", #Create a histogram of the sampling distribution
main="Approximate Sampling Distribution \n For Average Absolute r",
xlab = "Average Absolute r", ylab="Frequency",
xlim= range(min(samp.distr)-.01,AbsRObs+.01) )
abline (v=(Crit95r), col="red") #Plot the critical value as a line
points(AbsRObs,0, col="red", pch=19) #Plot the observed value point
hist(samp.distsig, freq=TRUE, col="cyan", #Create a histogram of the sampling distribution
main="Approximate Sampling Distribution \n For Number Significant",
xlab = "Number Statistically Significant", ylab="Frequency",
xlim= range(min(samp.distsig)-1,(SigObs+1)))
abline (v=(Crit95Sig), col="red") #Plot the critical value as a line
points(SigObs,0, col="red", pch=19) #Plot the observed value point
}
return(results)
}
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.