Nothing
######################################################################
#
# gtest.R
#
# copyright (c) 2004, Carter T. Butts <buttsc@uci.edu>
# Last Modified 10/5/20
# Licensed under the GNU General Public License version 2 (June, 1991)
# or later.
#
# Part of the R/sna package
#
# This file contains routines related to null hypothesis testing.
#
# Contents:
# cug.test
# cugtest
# plot.cug.test
# plot.cugtest
# plot.qaptest
# print.cug.test
# print.cugtest
# print.qaptest
# print.summary.cugtest
# print.summary.qaptest
# qaptest
# summary.cugtest
# summary.qaptest
#
######################################################################
#Second generation (univariate) CUG test routine
cug.test<-function(dat,FUN,mode=c("digraph","graph"),cmode=c("size","edges","dyad.census"),diag=FALSE,reps=1000,ignore.eval=TRUE,FUN.args=list()){
if(ignore.eval){ #Handle everything in edgelist form here
#Pre-process the raw input
dat<-as.edgelist.sna(dat)
if(is.list(dat))
return(sapply(dat,cug.test,FUN=FUN,mode=mode,cmode=cmode,diag=diag, reps=reps,ignore.eval=ignore.eval,FUN.args=FUN.args))
#End pre-processing
n<-attr(dat,"n")
if(!diag){
dat<-dat[dat[,1]!=dat[,2],,drop=FALSE]
attr(dat,"n")<-n
ndiag<-0
}else
ndiag<-sum(dat[,1]==dat[,2])
mode<-match.arg(mode)
cmode<-match.arg(cmode)
#Generate conditioning statistics
if(cmode=="size"){
m<-NULL
dc<-NULL
}else if(cmode=="edges"){
m<-switch(match.arg(mode),
graph=(NROW(dat)-ndiag)/2+ndiag,
digraph=NROW(dat)
)
dc<-NULL
}else if(cmode=="dyad.census"){
m<-NULL
dc<-dyad.census(dat)
}
#Generate randomization functions
getstat<-function(d){do.call(fun,c(list(d),FUN.args))}
drawrep<-switch(cmode,
size=function(n,...){rgraph(n,1,mode=mode,diag=diag,tprob=0.5, return.as.edgelist=TRUE)},
edges=function(n,m,...){rgnm(n=1,nv=n,m=m,mode=mode,diag=diag, return.as.edgelist=TRUE)},
dyad.census=function(n,dc,...){rguman(n=1,nv=n,mut=dc[1],asym=dc[2], null=dc[3],method="exact",return.as.edgelist=TRUE)},
)
}else{ #For valued data, we're going to use adjacency matrices (sorry)
#Pre-process the raw input
dat<-as.sociomatrix.sna(dat)
if(is.list(dat))
return(sapply(dat,cug.test,FUN=FUN,mode=mode,cmode=cmode,diag=diag, reps=reps,ignore.eval=ignore.eval,FUN.args=FUN.args))
else if(length(dim(dat))>2)
return(apply(dat,1,cug.test,FUN=FUN,mode=mode,cmode=cmode,diag=diag, reps=reps,ignore.eval=ignore.eval,FUN.args=FUN.args))
#End pre-processing
n<-NROW(dat)
m<-NULL
dc<-NULL
mode<-match.arg(mode)
cmode<-match.arg(cmode)
getstat<-function(d){do.call(fun,c(list(d),FUN.args))}
drawrep<-switch(cmode,
size=function(n,...){rgraph(n,1,mode=mode,diag=diag,tprob=0.5)},
edges=switch(mode,
digraph=function(n,...){g<-dat;
g[upper.tri(g,diag=diag)|lower.tri(g)]<- sample(g[upper.tri(g,diag=diag)|lower.tri(g)]); g
},
graph=function(n,...){g<-dat;
g[upper.tri(g,diag=diag)]<-sample(g[upper.tri(g,diag=diag)]);
g[lower.tri(g)]<-t(g)[lower.tri(g)]; g
}
),
dyad.census=function(n,...){g<-rewire.ud(dat,1)[1,,];
if(diag) diag(g)<-sample(diag(g)); g
},
)
}
#Set things up
fun<-match.fun(FUN)
if("mode"%in%names(formals(fun)))
callmode<-TRUE
else
callmode<-FALSE
if("diag"%in%names(formals(fun)))
calldiag<-TRUE
else
calldiag<-FALSE
if(callmode)
FUN.args$mode<-mode
if(calldiag)
FUN.args$diag<-diag
obs<-getstat(dat)
#Draw replicate stats
repstats<-vector()
for(i in 1:reps){
repstats[i]<-getstat(drawrep(n=n,m=m,dc=dc))
}
out<-list(obs.stat=obs,rep.stat=repstats,mode=mode,diag=diag,cmode=cmode, plteobs=mean(repstats<=obs),pgteobs=mean(repstats>=obs),reps=reps)
class(out)<-"cug.test"
out
}
#cugtest - Generate, print, and plot CUG (conditional uniform graph) test
#objects.
cugtest<-function(dat,FUN,reps=1000,gmode="digraph",cmode="density",diag=FALSE,g1=1,g2=2,...){
out<-list()
#Pre-process the raw input
dat<-as.sociomatrix.sna(dat)
#End pre-processing
#First, find the test value for fun on dat
fun<-match.fun(FUN)
out$testval<-fun(dat,g1=g1,g2=g2,...)
#Next, determine on what we are conditioning
if(cmode=="density"){
d<-c(gden(dat,g=g1,mode=gmode,diag=diag), gden(dat,g=g2,mode=gmode,diag=diag))
}else if(cmode=="ties"){
if(is.list(dat)){
tie1<-dat[[g1]]
tie2<-dat[[g2]]
}else{
tie1<-dat[g1,,]
tie2<-dat[g2,,]
}
}else{
d<-c(0.5,0.5)
}
if(is.list(dat)){ #Get the graph sizes
n1<-dim(dat[[g1]])[2]
n2<-dim(dat[[g2]])[2]
}else{
n1<-dim(dat)[2]
n2<-dim(dat)[2]
}
#Now, perform reps replications on random recreations of the data
out$dist<-vector(mode="numeric",length=reps)
for(i in 1:reps){
if(cmode=="ties"){ #Generate random replicates
dat1<-rgraph(n1,diag=diag,mode=gmode,tielist=tie1)
dat2<-rgraph(n2,diag=diag,mode=gmode,tielist=tie2)
}else{
dat1<-rgraph(n1,tprob=d[1],diag=diag,mode=gmode)
dat2<-rgraph(n2,tprob=d[2],diag=diag,mode=gmode)
}
if(n1==n2){ #Combine into single structure
datc<-array(dim=c(2,n1,n1))
datc[1,,]<-dat1
datc[2,,]<-dat2
}else
datc<-list(dat1,dat2)
out$dist[i]<-fun(datc,g1=1,g2=2,...) #Compute replicate stat
}
#Find p values
out$pgreq<-mean(as.numeric(out$dist>=out$testval))
out$pleeq<-mean(as.numeric(out$dist<=out$testval))
class(out)<-c("cugtest","cug")
out
}
#plot.cugtest - Plotting method for cugtest
plot.cug.test<-function(x,main="Univariate CUG Test", sub=paste("Conditioning:",x$cmode,"Reps:",x$reps),...){
xl<-range(c(x$rep.stat,x$obs.stat))
hist(x$rep.stat,xlim=xl,xlab="CUG Replicates",prob=TRUE,main=main,sub=sub,...)
abline(v=x$obs.stat,col=2,lwd=3)
}
#plot.cugtest - Plotting method for cugtest
plot.cugtest<-function(x,mode="density",...){
if(mode=="density"){
plot(density(x$dist),main="Estimated Density of CUG Replications",xlab="Test Statistic",...)
}else{
hist(x$dist,main="Histogram of CUG Replications",xlab="Test Statistic",...)
}
abline(v=x$testval,lty=2)
}
#plot.qaptest - Plotting method for qaptest
plot.qaptest<-function(x,mode="density",...){
if(mode=="density"){
plot(density(x$dist),main="Estimated Density of QAP Replications",xlab="Test Statistic",...)
}else{
hist(x$dist,main="Histogram of QAP Replications",xlab="Test Statistic",...)
}
abline(v=x$testval,lty=2)
}
#print.cug.test - Print method for cug.test
print.cug.test<-function(x,...){
cat("\nUnivariate Conditional Uniform Graph Test\n\n")
cat("Conditioning Method:",x$cmode,"\nGraph Type:",x$mode,"\nDiagonal Used:", x$diag,"\nReplications:",x$reps,"\n\n")
cat("Observed Value:",x$obs.stat,"\n")
cat("Pr(X>=Obs):",x$pgteobs,"\n")
cat("Pr(X<=Obs):",x$plteobs,"\n\n")
}
#print.cugtest - Print method for cugtest
print.cugtest<-function(x,...){
cat("\nCUG Test Results\n\n")
cat("Estimated p-values:\n")
cat("\tp(f(rnd) >= f(d)):",x$pgreq,"\n")
cat("\tp(f(rnd) <= f(d)):",x$pleeq,"\n\n")
}
#print.qaptest - Print method for qaptest
print.qaptest<-function(x,...){
cat("\nQAP Test Results\n\n")
cat("Estimated p-values:\n")
cat("\tp(f(perm) >= f(d)):",x$pgreq,"\n")
cat("\tp(f(perm) <= f(d)):",x$pleeq,"\n\n")
}
#print.summary.cugtest - Print method for summary.cugtest
print.summary.cugtest<-function(x,...){
cat("\nCUG Test Results\n\n")
cat("Estimated p-values:\n")
cat("\tp(f(rnd) >= f(d)):",x$pgreq,"\n")
cat("\tp(f(rnd) <= f(d)):",x$pleeq,"\n")
cat("\nTest Diagnostics:\n")
cat("\tTest Value (f(d)):",x$testval,"\n")
cat("\tReplications:",length(x$dist),"\n")
cat("\tDistribution Summary:\n")
cat("\t\tMin:\t",quantile(x$dist,probs=0,names=FALSE),"\n")
cat("\t\t1stQ:\t",quantile(x$dist,probs=0.25,names=FALSE),"\n")
cat("\t\tMed:\t",quantile(x$dist,probs=0.5,names=FALSE),"\n")
cat("\t\tMean:\t",mean(x$dist),"\n")
cat("\t\t3rdQ:\t",quantile(x$dist,probs=0.75,names=FALSE),"\n")
cat("\t\tMax:\t",quantile(x$dist,probs=1,names=FALSE),"\n")
cat("\n")
}
#print.summary.qaptest - Print method for summary.qaptest
print.summary.qaptest<-function(x,...){
cat("\nQAP Test Results\n\n")
cat("Estimated p-values:\n")
cat("\tp(f(perm) >= f(d)):",x$pgreq,"\n")
cat("\tp(f(perm) <= f(d)):",x$pleeq,"\n")
cat("\nTest Diagnostics:\n")
cat("\tTest Value (f(d)):",x$testval,"\n")
cat("\tReplications:",length(x$dist),"\n")
cat("\tDistribution Summary:\n")
cat("\t\tMin:\t",quantile(x$dist,probs=0,names=FALSE),"\n")
cat("\t\t1stQ:\t",quantile(x$dist,probs=0.25,names=FALSE),"\n")
cat("\t\tMed:\t",quantile(x$dist,probs=0.5,names=FALSE),"\n")
cat("\t\tMean:\t",mean(x$dist),"\n")
cat("\t\t3rdQ:\t",quantile(x$dist,probs=0.75,names=FALSE),"\n")
cat("\t\tMax:\t",quantile(x$dist,probs=1,names=FALSE),"\n")
cat("\n")
}
#qaptest - Generate a QAP test object
qaptest<-function(dat,FUN,reps=1000,...){
out<-list()
#Coerce to sociomatrix form - rmperm is going to do that anyway....
dat<-as.sociomatrix(dat)
#First, find the test value for fun on dat
fun<-match.fun(FUN)
out$testval<-fun(dat,...)
#Now, perform reps replications on random permutations of the data
out$dist<-vector(mode="numeric",length=reps)
for(i in 1:reps){
out$dist[i]<-fun(rmperm(dat),...)
}
#Find p values
out$pgreq<-mean(as.numeric(out$dist>=out$testval))
out$pleeq<-mean(as.numeric(out$dist<=out$testval))
class(out)<-c("qaptest","qap")
out
}
#summary.cugtest - Summary method for cugtest
summary.cugtest<-function(object, ...){
out<-object
class(out)<-c("summary.cugtest",class(out))
out
}
#summary.qaptest - Summary method for qaptest
summary.qaptest<-function(object, ...){
out<-object
class(out)<-c("summary.qaptest",class(out))
out
}
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.