############################### Study2Template.txt ################################### # # # Study2 is on F-tests 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(86). 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. Say this is a one-way ANOVA, so # p = df1+1, and df2 = n-p = n-df1-1 geq 1, so n > df1+2 # This function is different in many places for F and chi-squared fun = function(es,wantpow,dfree1,lambda=86) { alpha = 0.05; # Sample size AFTER selection is Poisson top = round(lambda + 6*sqrt(lambda)) n = (dfree1+2):top # Sum Poisson up to top (6 sigma above mean) dfree2 = n-dfree1-1 cv=qf(1-alpha,dfree1,dfree2) epow = sum((1-pf(cv,df1=dfree1,df2=dfree2,ncp=n*es^2))*dpois(n,lambda)) # F # 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,dfree1=1,wantpow=0.25)$par EffectSize[1,3] = nlminb(start=0.01, objective=fun,lower=0,dfree1=1,wantpow=0.50)$par EffectSize[1,4] = nlminb(start=.5, objective=fun,lower=0,dfree1=1,wantpow=0.75)$par EffectSize[2,2] = nlminb(start=0.01, objective=fun,lower=0,dfree1=3,wantpow=0.25)$par EffectSize[2,3] = nlminb(start=0.01, objective=fun,lower=0,dfree1=3,wantpow=0.50)$par EffectSize[2,4] = nlminb(start=0.6, objective=fun,lower=0,dfree1=3,wantpow=0.75)$par EffectSize[3,2] = nlminb(start=0.01, objective=fun,lower=0,dfree1=5,wantpow=0.25)$par EffectSize[3,3] = nlminb(start=0.01, objective=fun,lower=0,dfree1=5,wantpow=0.50)$par EffectSize[3,4] = nlminb(start=0.6, objective=fun,lower=0,dfree1=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. DFree1 = c(1,3,5); index = 1; Dfree1 = DFree1[index] # Lambda = Poisson mean sample size was set above. # cv needs to be set in the loop -- varies with N nsim=25 # 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) DF2mat = Nmat-Dfree1-1 crit = qf(0.95,df1=Dfree1,df2=DF2mat) NCPmat = es^2*Nmat Powmat = 1-pf(crit,df1=Dfree1,df2=DF2mat,ncp=NCPmat) Ymat = rsigF(nsim*k,DF1=Dfree1,DF2=DF2mat,NCP=NCPmat) # Format Ymat, Nmat Powmat to process by row # Format Nmat and Ymat to process by row dim(Ymat) = dim(Powmat) = dim(Nmat) = dim(DF2mat) = c(nsim,k) for(Simulation in 1:nsim) { N = Nmat[Simulation,]; Y = Ymat[Simulation,]; DF2 = DF2mat[Simulation,] TruePower = Powmat[Simulation,] # k values, all potentially different MeanTruePower = mean(TruePower) # Sample mean of True Power values rowno = rowno+1 Pcurve = heteroNpcurveF(FF=Y,dfree1=Dfree1,dfree2=DF2) Puniform = heteroNpunifF(FF=Y,dfree1=Dfree1,dfree2=DF2,CI=F) MaxLike = heteroNmleF(FF=Y,dfree1=Dfree1,dfree2=DF2,CI=F) pval = 1-pf(Y,df1=Dfree1,df2=DF2) Zcurve = zcurve(pval,Plot=0,Verbose=F,ncomp=5) resultz = c(Pcurve, Puniform, MaxLike, Zcurve) # Put them in the table simdata[rowno,] = c(Simulation,Nstudies[i], AvePower[j], MeanTruePower, resultz) } # 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) # But to accomodate large data files, ... # 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?