# demo/powersim.R In highD2pop: Two-Sample Tests for Equality of Means in High Dimension

```##### Set simulation parameters here.

simparms <- list	(
n = c(30,45),
m = c(45,60),
p = 100,
beta = c(0,.1,.2,.4,.6,.8,.9,1),
delta = .5,
muX = function(p) return(rep(0,p)),
muY = function(p,b,delta) {
times1 <- b * p
times2 <- (1-b)*p+1e-7
mu <- c(rep(delta,times1),rep(0,times2))
return(mu)
},
commoncov = TRUE,
VarScaleY = 2,
dep = c('IND','ARMA','LR'),
ARMAparms = list(coefs=list(ma=c(.2,.3) , ar=c(.4,-.1))),
LRparm = .75,
S = 25,
#innov = function(n,...) rdblepareto(n,16.5,8),
#innov = function(n,...) rnorm(n,0,1),
innov = function(n,...) rgammashift(n,4,2),
heteroscedastic=FALSE,
het.diag = function(p) return(diag(1/2 + rexp(p,2)))
)

testparms <- list(	r = 15,
smoother = c("parzen","trapezoid"),
ntoorderminus=c(0,1,2))

library(highD2pop)

highD2popBetaPower.sim <- function 	(
ChenQin = FALSE,
SK	= FALSE,
CLX = FALSE,
GCT = TRUE,
simparms,
testparms
)
{
# Hmmm, I may just have to redefine the function if I change my simparms...

output<-list
output<-list(output,output,output,output,output,output,output,output)  # l
output<-list(output,output,output)	# k
output<-list(output)	# j
output<-list(output,output) 	# i

for(i in 1:length(simparms\$n)){
for(j in 1:length(simparms\$p)){
for(k in 1:length(simparms\$dep)){
for(l in 1:length(simparms\$beta)){

i<-1
j<-1
k<-1
l<-3

DATA <-	build2popData(
n = simparms\$n[i],
m = simparms\$m[i],
p = simparms\$p,
muX = simparms\$muX(simparms\$p),
muY = simparms\$muY(simparms\$p,simparms\$beta[l],simparms\$delta),
commoncov = simparms\$commoncov,
VarScaleY = simparms\$VarScaleY,
dep = simparms\$dep[k],
ARMAparms = simparms\$ARMAparms,
LRparm = simparms\$LRparm,
S = simparms\$S,
innov = simparms\$innov,
heteroscedastic=simparms\$heteroscedastic,
het.diag = simparms\$het.diag(simparms\$p)
)

if(SK==TRUE) SK.results <- SK.sim(DATA) else SK.results=NULL
if(ChenQin==TRUE) ChenQin.results <- ChenQin.sim(DATA) else ChenQin.results=NULL
if(CLX==TRUE) CLX.results <- CLX.sim.Covtest(DATA) else CLX.results=NULL
if(GCT==TRUE) GCT.results11 <- GCT.sim(DATA,r = testparms\$r[1],smoother=testparms\$smoother[1],ntoorderminus=testparms\$ntoorderminus[1]) else GCT.results1=NULL
if(GCT==TRUE) GCT.results12 <- GCT.sim(DATA,r = testparms\$r[1],smoother=testparms\$smoother[1],ntoorderminus=testparms\$ntoorderminus[2]) else GCT.results2=NULL
if(GCT==TRUE) GCT.results13 <- GCT.sim(DATA,r = testparms\$r[1],smoother=testparms\$smoother[1],ntoorderminus=testparms\$ntoorderminus[3]) else GCT.results3=NULL
if(GCT==TRUE) GCT.results21 <- GCT.sim(DATA,r = testparms\$r[1],smoother=testparms\$smoother[2],ntoorderminus=testparms\$ntoorderminus[1]) else GCT.results1=NULL
if(GCT==TRUE) GCT.results22 <- GCT.sim(DATA,r = testparms\$r[1],smoother=testparms\$smoother[2],ntoorderminus=testparms\$ntoorderminus[2]) else GCT.results2=NULL
if(GCT==TRUE) GCT.results23 <- GCT.sim(DATA,r = testparms\$r[1],smoother=testparms\$smoother[2],ntoorderminus=testparms\$ntoorderminus[3]) else GCT.results3=NULL

output[[i]][[j]][[k]][[l]] <- list(	SK.results = SK.results,
ChenQin.results = ChenQin.results,
CLX.results = CLX.results,
GCT.results11 = GCT.results11,
GCT.results12 = GCT.results12,
GCT.results13 = GCT.results13,
GCT.results21 = GCT.results21,
GCT.results22 = GCT.results22,
GCT.results23 = GCT.results23,
n = simparms\$n[i],
m = simparms\$m[i],
p = simparms\$p,
muX = simparms\$muX(simparms\$p),
muY = simparms\$muY(simparms\$p,simparms\$beta[l],simparms\$delta),
commoncov = simparms\$commoncov,
VarScaleY = simparms\$VarScaleY,
dep = simparms\$dep[k],
ARMAparms = simparms\$ARMAparms,
LRparm = simparms\$LRparm,
S = simparms\$S,
innov = simparms\$innov,
heteroscedastic=simparms\$heteroscedastic,
het.diag = simparms\$het.diag(simparms\$p)
)
print(paste(i,j,k,l))

}
}
}
}

return(output)

}

output <- highD2popBetaPower.sim(SK=TRUE,ChenQin=TRUE,CLX=TRUE,GCT=TRUE,simparms=simparms,testparms=testparms)

##### Plot the results.

par(mfrow = c(2,3))

par(oma = c(2.5,2.5,4,1),
mar = c(0,0,0,0),
lwd = 1,
cex.axis = .8,
cex.lab = .85,
tck = -.02)

Pzen0 <- Pzen2 <- ChQ <- SK <- CLX <- matrix(0,length(output),length(output[[1]][[1]][[1]]))

for( i in 1:length(output))
for(j in 1:3)
{
for(l in 1:length(output[[1]][[1]][[1]]))
{
Pzen0[i,l] <- mean(output[[i]][[1]][[j]][[l]]\$GCT.results11[,2] < .05) # parzen, r = 20
Pzen2[i,l] <- mean(output[[i]][[1]][[j]][[l]]\$GCT.results13[,2] < .05) # parzen, r = 20
ChQ[i,l] <- mean(output[[i]][[1]][[j]][[l]]\$ChenQin.results[,2] < .05) # Ch-Q
SK[i,l] <- mean(output[[i]][[1]][[j]][[l]]\$SK.results[,2] < .05) # SK
CLX[i,l] <- mean(output[[i]][[1]][[j]][[l]]\$CLX.results[,2] < .05) # CLX
}

plot(	Pzen0[i,],
ylim = c(0,1),
xlim = c(1,8),
yaxt = "n",
ylab = "",
type = "l",
xaxt = "n",
xlab = "",
lwd = 1)

lines(Pzen0[i,],lwd=1)
lines(Pzen2[i,],lty=5,lwd=1)
lines(ChQ[i,],lty=3,lwd=1)
lines(SK[i,],lty=2,lwd=1)
lines(CLX[i,],lty=4,lwd=1)

if(i==1 & j==2){

x.pos <- grconvertX(.5,from="nfc",to="user")
y.pos <- grconvertY(1.125,from="nfc",to="user")

legend(	x=x.pos,
y=y.pos,
lwd=c(1,1,1,1,1),
lty=c(1,5,3,2,4),
legend=c("mod-p GCT","lg-p GCT","Ch-Q","SK","CLX"),
bty="n",
seg.len=3,
horiz = TRUE,
xpd=NA,
xjust = .5,
cex=1.25)
}

if(i==1)
{

text(x=1,y=.5,"(n,m) = (45,60)",srt=90,cex=.85)

}

if(i==2)
{
text(x=1,y=.5,"(n,m) = (90,120)",srt=90,cex=.85)
}

abline(h=.05,lwd=.5)
text(x=7.5,y=.025, expression(alpha==.05),cex=.6)

}

text(x=7.5,y=.025, expression(alpha==.05),cex=.6)

mtext(outer=TRUE,side=1,expression(beta), line=1.5,cex=1)
mtext(outer=TRUE,side=2,"Power",line=1.25,cex=1)
mtext(outer=TRUE,side=3,text="Power curves under skewed innovations and equal covariance matrices",line=2.5,cex=1)
```

## Try the highD2pop package in your browser

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

highD2pop documentation built on May 29, 2017, 6:09 p.m.