############################### Study1Template.txt ################################### # # # Study1 is a chi-squared test with significance level 0.05 under conditions # # of publication bias and heterogeneity in sample size only. In this version, # # sample size after selection is Poisson(80). The 7 variables are # # Simulation Nstudies AvePower MeanTruePower Pcurve, Puniform, MaxLike, Zcurve # # AvePower is population mean power, while MeanTruePower is SAMPLE mean true power. # # Copy to template.txt, and run with runR on Depth Charge or sniffR on PsychLab # # # # # # Note the degree of freedom could be 1, 3 or 5 in separate runs. # # # ###################################################################################### rm(list=ls()) options(scipen=999) # To avoid scientific notation source("http://www.utstat.toronto.edu/~brunner/Rfunctions/estimatR.txt") Lambda = 86 # Poisson mean sample size # PsychScience = scan("http://www.utstat.toronto.edu/~brunner/data/power/PsychScience.urn.txt") # # Just n of 20 or more # set.seed(888) # PsychScience = sample(subset(PsychScience,PsychScience>19)); # length(PsychScience) # 6444 # PsychScience = PsychScience[1:6000]; mean(PsychScience) # 86.04383 # rm(.Random.seed) ############### Find needed Effect sizes ############### # Minimize this to find effect size giving desired power # after selection for significance fun = function(es,wantpow,dfreedom,lambda=Lambda) { alpha = 0.05; cv=qchisq(1-alpha,dfreedom) # Sample size AFTER selection is Poisson top = round(Lambda + 6*sqrt(Lambda)) n = 1:top # Sum Poisson up to top (6 sigma above mean) epow = sum((1-pchisq(cv,df=dfreedom,ncp=n*es^2))*dpois(n,Lambda)) # cat("es = ",es," Expected power = ",epow,"\n") (epow-wantpow)^2 } # End of all the fun EffectSize = matrix(0,3,4) rownames(EffectSize) = c(1,3,5); colnames(EffectSize) = c(0,0.25,0.50,0.75) # EffectSize EffectSize[1,2] = nlminb(start=0.01, objective=fun,lower=0,df=1,wantpow=0.25)$par EffectSize[1,3] = nlminb(start=0.01, objective=fun,lower=0,df=1,wantpow=0.50)$par EffectSize[1,4] = nlminb(start=.5, objective=fun,lower=0,df=1,wantpow=0.75)$par EffectSize[2,2] = nlminb(start=0.01, objective=fun,lower=0,df=3,wantpow=0.25)$par EffectSize[2,3] = nlminb(start=0.01, objective=fun,lower=0,df=3,wantpow=0.50)$par EffectSize[2,4] = nlminb(start=0.6, objective=fun,lower=0,df=3,wantpow=0.75)$par EffectSize[3,2] = nlminb(start=0.01, objective=fun,lower=0,df=5,wantpow=0.25)$par EffectSize[3,3] = nlminb(start=0.01, objective=fun,lower=0,df=5,wantpow=0.50)$par EffectSize[3,4] = nlminb(start=0.6, objective=fun,lower=0,df=5,wantpow=0.75)$par cat("\n Effect sizes needed for desired population mean power \n") cat(" Rows are df, columns are power values \n\n") print(EffectSize) ####################### # Setup # Run separately for each df. Change name of output file accordingly. DFREE = c(1,3,5); index = 1; Dfree = DFREE[index] # Lambda = Poisson mean sample size was set above. cv = qchisq(0.95,Dfree) nsim=1250 # Number of simulated data sets in each condition: 1,250 * 8 = 10,000 siglevel = 0.05 # Loop through these variables Nstudies = c(15,25,50,100,250) # 5 levels AvePower = c(0.05,0.25,0.50,0.75) # 4 levels nrows = 5*4*nsim; simdata = matrix(0,nrows,8) # Simulation number plus seven variables Methods = c("Pcurve","Puniform","MaxLike","Zcurve") # Make the main simulation a function so it's easier to time when run from the unix # command line. doit = function() { rowno = 0 # Getting ready to replace one row at a time, much faster than rbind. for(i in 1:length(Nstudies)) { k = Nstudies[i] for(j in 1:length(AvePower)) { es = EffectSize[index,j] # Generate data all at once in matrices to process one row at a time. Nmat = rpois(nsim*k,Lambda) NCPmat = es^2*Nmat; Powmat = 1-pchisq(cv,df=Dfree,ncp=NCPmat) Ymat = rsigCHI(nsim*k,DF=Dfree,NCP=NCPmat) # Format Ymat, Nmat Powmat to process by row dim(Ymat) = dim(Powmat) = dim(Nmat) = c(nsim,k) # Format Nmat and Ymat to process by row dim(Nmat) = c(nsim,k); dim(Ymat) = c(nsim,k); dim(Powmat) = c(nsim,k) for(Simulation in 1:nsim) { N = Nmat[Simulation,]; Y = Ymat[Simulation,] TruePower = Powmat[Simulation,] # k values, all potentially different MeanTruePower = mean(TruePower) # Sample mean of True Power values rowno = rowno+1 Pcurve = heteroNpcurveCHI(Y=Y,dfree=Dfree,nn=N) Puniform = heteroNpunifCHI(Y=Y,dfree=Dfree,nn=N,CI=F) MaxLike = heteroNmleCHI(Y=Y,dfree=Dfree,nn=N,CI=F) pval = 1-pchisq(Y,df=Dfree) Zcurve = zcurve(pval,Plot=0,Verbose=F) simdata[rowno,] = c(Simulation,Nstudies[i], AvePower[j], MeanTruePower, Pcurve, Puniform, MaxLike, Zcurve) } # Next Simulation } # Next value of TruePower } # Next value of Nstudies # Write the data frame to a file colnames(simdata) = c("Simulation", "Nstudies", "PopMeanPower", "SampleMeanPower", "Pcurve", "Puniform", "MaxLike", "Zcurve") # Write the data frame to a file. File will not have names. Could do # write.table(data.frame(simdata),"Study1.1.Data.txt",quote=F) # Data are written one column at a time, so write the transpose write(t(simdata),"simoutX.txt",ncolumns=8) # The only appearance of upper case x. } # End function doit, which just does it. system.time(doit()) # How long did it take to do it?