If(exists('compare - Horvath Lab UCLA



Supplement to Network Edge Orienting:

Study of Single Anchor vs. Orthogonal Causal Anchor LEO scoring: statistical power and false positive rate graphs

Jason E. Aten

5 October 2007

# here is the archival versions of the two primary plots we are

# producing here: 2 encapsulated postscript files inside this zip file.

[pic]

# double click to unzip these postscript graphs and

# save them into your manuscript directory.

# Now to reproduce the generation of these figures:

# Set the working directory.

setwd("c:/good")

# Source network functions

source("neo.txt")

# unzip this data file into your working directory.

[pic]

# Function to grid sample a range of parameter values to be simulated on

# a large scale. This was used to prepare to simulate the data contained

# in scores424.rdat.

if(exists("generate.simulation.params") ) rm(generate.simulation.params);

generate.simulation.params=function() {

min.hb = min.ha = .2

max.hb = max.ha = .6

# min.hb = min.ha = .5

# max.hb = max.ha = .5

num.h.points= 3

range.ha = max.ha - min.ha

range.hb = max.hb - min.hb

n.sample.points = 3

omega.n.sample.points = 6

min.gamma2=0

res=data.frame()

ha.choices = seq(min.ha,max.ha,range.ha/num.h.points);

for (ha in ha.choices) {

#ha=ha.choices[1]

hb.choices = seq(min.hb,max.hb,range.hb/num.h.points);

for (hb in hb.choices) {

#hb=hb.choices[1]

min.errB = .05

min.errA = .05

= .7

omega.max = round(min(sqrt(1-ha-min.errA),),3)

omega.min = 0

omega.range = omega.max-omega.min;

omega.choices = seq(omega.min, omega.max, omega.range/omega.n.sample.points)

for (omega in omega.choices) {

percent.max.omega = round(omega/omega.max,3)

max.gamma2 = round(min((1-ha-omega^2-min.errA)/(1+2*omega),1-hb-min.errB,.6),3)

# max.gamma2 = round(min(1-hb-min.errB, 1-ha- min.errA, .6),3)

gamma2.choices = round(seq(min.gamma2,max.gamma2,(max.gamma2-min.gamma2)/n.sample.points),3);

for (gamma2 in gamma2.choices) {

# gamma2=gamma2.choices[1]

gamma=sqrt(gamma2)

eA=1-2*omega*gamma^2 - gamma^2 - omega^2 - ha;

if (eA < 0) { stop("algebra error.") }

eB=1-gamma2 - hb;

if (eB < 0) { stop("algebra error.") }

rho.ab=gamma2+omega

if (rho.ab>1) { stop("algebra error.") }

q=data.frame(ha,hb,gamma2,omega,eA,eB,rho.ab,percent.max.omega)

res=rbind(res,q)

}

}

}

}

round(res,3)

}

# pull out the model probability

if(exists("extract.chisq.right.log10prob") ) rm(extract.chisq.right.log10prob);

extract.chisq.right.log10prob=function(sem.obj) {

if (!inherits(sem.obj,"summary.sem")) stop("bad object passed to extract.chisq.right.log10prob");

-pchisq(sem.obj$chisq,sem.obj$df,lower.tail=FALSE,log.p=TRUE)*log.to.log10

}

# ===================================================

# this function computes the standard error

if (exists("stderr1")) rm(stderr1)

stderr1 th1,1,0)

LEO.Icut2=ifelse(score2> th2,1,0)

# group data by percent.max.omega, then get the frequency of rmsea being good in each group

# ; the mean or expected value of an indicator is just the probability or frequency of that

# indicator.

# tapply automatically sorts the percent.max.omega groups into ascending order.

# for 1

POWER1=as.vector(tapply(LEO.Icut1,scores424$percent.max.omega[w1],mean));

group.names1=round(sort(unique(scores424$percent.max.omega[w1])),2)

names(POWER1)=as.character(group.names1)

SE.POWER1 = as.vector(tapply(LEO.Icut1,scores424$percent.max.omega[w1],stderr1))

#ci.l1 = POWER1-1.96*SE.POWER1

#ci.u1 = POWER1+1.96*SE.POWER1

ci.l1 = POWER1-SE.POWER1

ci.u1 = POWER1+SE.POWER1

# for 2

POWER2=as.vector(tapply(LEO.Icut2,scores424$percent.max.omega[w2],mean));

group.names2=round(sort(unique(scores424$percent.max.omega[w2])),2)

names(POWER2)=as.character(group.names2)

SE.POWER2 = as.vector(tapply(LEO.Icut2,scores424$percent.max.omega[w2],stderr1))

#ci.l2 = POWER2-1.96*SE.POWER2

#ci.u2 = POWER2+1.96*SE.POWER2

ci.l2 = POWER2-SE.POWER2

ci.u2 = POWER2+SE.POWER2

par(mfrow=c(1,1))

if (!do.postscript) windows()

if (do.postscript) postscript(file=paste(sep="",gsub(" ",".",score1.name),".vs.",gsub(" ",".",score2.name),".BARPLOT.vs.percent.true.signal.side.by.side.",hahb1,".vs.",hahb2,".ps"),horizontal=FALSE)

mp1 ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download