daily pastebin goal
47%
SHARE
TWEET

R simulation

a guest Jan 30th, 2013 38 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. require(fields) #Used for image.plot
  2.  
  3.  
  4. ###Detect RStudio (modifies plot settings:no live plot)
  5. using.RStudio=nzchar(Sys.getenv("RSTUDIO_USER_IDENTITY"))
  6.  
  7.  
  8.  
  9. ###General Settings###
  10. runs<-1000                      # Number of Simulations to Run
  11. max.samps<-15                   # Maximum number of Samples to test
  12. initial.samp.size<-5            # must be at least 2
  13. samp.interval<-5                        # must be at least 1
  14. effect<-0  # Difference between population means
  15. pvalue<-.05
  16.  
  17. #Live Plot Settings
  18. live.plot=F                     # If false, p values traces from last 10 simulations will be plotted
  19. sleep.time<-.1          # In seconds, use to slow down live plot
  20. live.plot.in.own.window=F       # If live.plot=T, also plot p value trace in its own device
  21.  
  22.  
  23. ###Investigator Bias Settings###
  24. ##Sequential Hypothesis Testing
  25. # Set true (T) to add samples to the previous sample.
  26. # Otherwise (set to F) each sample is independant
  27. sequential=T           
  28.  
  29. ##Multiple Statistical Tests
  30. # If True, choose lowest of t-test and Mann-Whitney U-test p-values.
  31. # Otherwise only perform t-test
  32. also.perform.Utest=T   
  33.  
  34. ##Multiple Outcomes/Publication Bias, e.g. only report significant comparisons
  35. #Set greater than 1 if multiple outcome measures will be tested.
  36. #Probabilities and effect sizes will be for any of the tests significant at each sample size
  37. #History of each is saved for sequential testing
  38. number.of.outcomes= 3
  39.  
  40. ##Drop Outliers
  41. # Set to true to drop outliers as defined by outlier.cutoff
  42. # In case of sequential, outliers are kept for next iteration
  43. # Outlier bias are multipliers of outlier cutoff for treatment group
  44. # Set > 1 to be less strict
  45. # Set < 1 for more strict
  46. drop.outliers=T
  47. outlier.method="MAD" # Defaults to "MAD" (median absolute deviation), can also use "SD" (number of standard deviations from mean)
  48. outlier.cutoff<-5               # Number of sample deviations from sample mean/median (upper and lower)
  49. outlier.bias.high<-1    # Bias multiplier (see above)
  50. outlier.bias.low<-1   # Bias multiplier (see above)
  51. n.cutoff<-2   # Replace values after dropping if sample size is less than the cutoff, min=2, max=sample size
  52.  
  53.  
  54. ###Population Distributions###
  55. #Treatment Group
  56. treatment.sd1<-1                # Treatment Group Primary Standard Deviation
  57. treatment.shape1<-2     # Treatment Group Primary Shape Parameter (set to 2 for normal)
  58. treatment.sd2<-1                # Treatment Group Secondary Standard Deviation (for bimodal)
  59. treatment.shape2<-2     # Treatment Group Secondary Shape Parameter (set to 2 for normal, for bimodal)
  60.  
  61. #Control Group
  62. control.sd1<-1          # Control Group Primary Standard Deviation
  63. control.shape1<-2               # Control Group Primary Shape Parameter (set to 2 for normal)
  64. control.sd2<-1          # Control Group Secondary Standard Deviation (for bimodal)
  65. control.shape2<-2               # Control Group Secondary Shape Parameter (set to 2 for normal, for bimodal)
  66.  
  67.  
  68.  
  69. ###Bimodal Settings###
  70. larger.bimodal.proportion<- .5  #Set to 1 for unimodal
  71. bimodal.offset<-8               #Difference Between Primary and Secondary Means for bimodal Populations
  72. bimodal.treatment.only=F        # Simulate treatment effect on subset of treatment group, control will be unimodal.
  73.  
  74.  
  75.  
  76.  
  77. ### Advanced: Random Parameter Generation Settings ###
  78. random.pop.parameters=F         # Generate Population Parameters (ignores settings above)
  79. random.bias=F                   # Generate Bias Parameters (ignores settings above)
  80. equal.pop.parameters=T          # Sets Treatment and Control Group Population Parameters equal
  81. random.shape.only=T             # Use sd settings above
  82. batch.size=1                    # Set >1 to run batch of random simulations
  83.  
  84.  
  85.  
  86.  
  87.  
  88. ###About Shape (p) Parameter###
  89. #Generalized normal distribution with Shape(p) Parameter#
  90. #Modified from package "normalp" to return numeric(0) rather than errors
  91. rnormp2<-function (n, mu = 0, sigmap = 1, p = 2, method = c("def", "chiodi"))
  92. {
  93.   if (!is.numeric(n) || !is.numeric(mu) || !is.numeric(sigmap) ||
  94.     !is.numeric(p))
  95.     return(numeric(0))
  96.   if (p < 1)
  97.     stop("p must be at least equal to one")
  98.   if (sigmap <= 0)
  99.     return(numeric(0))
  100.   method <- match.arg(method)
  101.   if (method == "def") {
  102.     qg <- rgamma(n, shape = 1/p, scale = p)
  103.     z <- qg^(1/p)
  104.     z <- ifelse(runif(n) < 0.5, -z, z)
  105.     x <- mu + z * sigmap
  106.   }
  107.   if (method == "chiodi") {
  108.     i <- 0
  109.     x <- c(rep(0, n))
  110.     while (i < n) {
  111.       u <- runif(1, -1, 1)
  112.       v <- runif(1, -1, 1)
  113.       z <- abs(u)^p + abs(v)^(p/(p - 1))
  114.       if (z <= 1) {
  115.         i <- i + 1
  116.         x[i] <- u * (-p * log(z)/z)^(1/p)
  117.         x[i] <- mu + x[i] * sigmap
  118.       }
  119.     }
  120.   }
  121.   x
  122. }
  123. # shape < 1 : invalid #
  124. #plot(seq(-5,5,by=.01),dnormp(seq(-5,5,by=.01), mu=0, sigmap=1, p=.5))
  125. #lines(seq(-5,5,by=.01), dnorm(seq(-5,5,by=.01), 0,1), lwd=2, col="Blue")
  126.  
  127. # 1<= shape < 2 : heavier tails, 1 is laplace dist #
  128. #plot(seq(-5,5,by=.01),dnormp(seq(-5,5,by=.01), mu=0, sigmap=1, p=1.5))
  129. #lines(seq(-5,5,by=.01), dnorm(seq(-5,5,by=.01), 0,1), lwd=2, col="Blue")
  130.  
  131. # shape = 2 : normal dist #
  132. #plot(seq(-5,5,by=.01),dnormp(seq(-5,5,by=.01), mu=0, sigmap=1, p=2))
  133. #lines(seq(-5,5,by=.01), dnorm(seq(-5,5,by=.01), 0,1), lwd=2, col="Blue")
  134.  
  135. # shape > 2 :lighter tails, inf is uniform dist #
  136. #plot(seq(-5,5,by=.01),dnormp(seq(-5,5,by=.01), mu=0, sigmap=1, p=10))
  137. #lines(seq(-5,5,by=.01), dnorm(seq(-5,5,by=.01), 0,1), lwd=2, col="Blue")
  138.  
  139.  
  140.  
  141.  
  142.  
  143. ###Start Script###
  144.  
  145. if (batch.size > 1){
  146.   batchmode=T
  147.   batch.results=NULL
  148.   batch.params=NULL
  149. }else{
  150.   batchmode=F
  151. }
  152.  
  153. for(b in 1:batch.size){
  154.  
  155.  
  156.   ###Randomly Set Bias###
  157.  
  158.   if(random.bias==T){
  159.     ###Investigator Bias Settings###
  160.     ##Sequential Hypothesis Testing
  161.     seq.samp<-round(runif(1,0,1))
  162.     if(seq.samp==1){
  163.       sequential=T
  164.     }else{
  165.       sequential=F
  166.     }          
  167.    
  168.     ##Multiple Statistical Tests
  169.     Utest.samp<-round(runif(1,0,1))
  170.     if(Utest.samp==1){
  171.       also.perform.Utest=T     
  172.     }else{
  173.       also.perform.Utest=F
  174.     }
  175.    
  176.     ##Multiple Outcomes/Publication Bias, e.g. only report significant comparisons
  177.     number.of.outcomes= round(runif(1, 1, 3),0)
  178.    
  179.     ##Drop Outliers
  180.     outlier.samp<-round(runif(1,0,1))
  181.     if(outlier.samp==1){
  182.       drop.outliers=T
  183.     }else{
  184.       drop.outliers=F
  185.     }
  186.     outlier.cutoff<-runif(1.5,1000)    
  187.     outlier.bias.high<-runif(1,.5,1.5) 
  188.     outlier.bias.low<-runif(1,.5,1.5)
  189.    
  190.    
  191.   } #End Random Bias Settings
  192.  
  193.  
  194.   ###Randomly Set Population Parameters###
  195.   if(random.pop.parameters==T){
  196.    
  197.     ###Population Distributions###
  198.     #Treatment Group
  199.     if(random.shape.only!=T){
  200.       treatment.sd1<-runif(1, .1, 5)
  201.     }          
  202.     treatment.shape1<-runif(1, 1, 10)
  203.     if(random.shape.only!=T){  
  204.       treatment.sd2<-runif(1, .1, 5)
  205.     }          
  206.     treatment.shape2<-runif(1, 1, 10)  
  207.    
  208.     #Control Group
  209.     if(equal.pop.parameters==T){
  210.       if(random.shape.only!=T){
  211.         control.sd1<-treatment.sd1
  212.       }        
  213.       control.shape1<-treatment.shape1
  214.       if(random.shape.only!=T){
  215.         control.sd2<-treatment.sd2
  216.       }        
  217.       control.shape2<-treatment.shape2 
  218.     }else{
  219.       if(random.shape.only!=T){
  220.         control.sd1<-runif(1, .1, 5)
  221.       }        
  222.       control.shape1<-runif(1, 1, 10)
  223.       if(random.shape.only!=T){
  224.         control.sd2<-runif(1, .1, 5)
  225.       }        
  226.       control.shape2<-runif(1, 1, 10)          
  227.     }
  228.    
  229.     ###Bimodal Settings###
  230.     larger.bimodal.proportion<- runif(1,.5,1)  
  231.     bimodal.offset<-runif(1, .1, 10)           
  232.    
  233.     if(equal.pop.parameters==T){
  234.       bimodal.treatment.only=F
  235.     }else{
  236.      
  237.       bimod.samp<-round(runif(1,0,1))
  238.      
  239.       if(bimod.samp==1){       
  240.         bimodal.treatment.only=T
  241.       }else{
  242.         bimodal.treatment.only=F
  243.       }
  244.      
  245.     }  
  246.    
  247.   }     #End Random Population Parameter Generation
  248.  
  249.  
  250.   ###Start Sampling###
  251.   if(initial.samp.size<2|| samp.interval<1){
  252.     print("Simulation not Run: minimum sample size =2, minimum sample interval=1")
  253.   }else{
  254.    
  255.    
  256.     ###Create Population Histogram and Live Chart###
  257.     smaller.bimodal.proportion<- (1-larger.bimodal.proportion)
  258.     bimodal.pop<-c(rep(1,larger.bimodal.proportion*1000),
  259.                    rep(0,smaller.bimodal.proportion*1000))
  260.    
  261.     if(larger.bimodal.proportion==1){
  262.       treatment.sd2<-NA
  263.       treatment.shape2<-NA
  264.       control.sd2<-NA
  265.       control.shape2<-NA
  266.       bimodal.offset<-NA
  267.       bimodal.treatment.only=F
  268.     }
  269.    
  270.     if(drop.outliers!=T){
  271.       outlier.cutoff<-NA
  272.       outlier.bias.high<-NA  
  273.       outlier.bias.low<-NA
  274.     }
  275.    
  276.     if(outlier.method!="MAD" && outlier.method!="SD"){
  277.       outlier.method="MAD"
  278.     }
  279.    
  280.     if(bimodal.treatment.only==T){
  281.       control.sd2<-NA
  282.     }
  283.    
  284.     if(bimodal.treatment.only==T){
  285.       control<-rnormp2(100000,0,control.sd1,control.shape1)
  286.     }else{
  287.       control<-c(rnormp2(larger.bimodal.proportion*100000,0,control.sd1,control.shape1),
  288.                  rnormp2(smaller.bimodal.proportion*100000,0+bimodal.offset,control.sd2,control.shape2))
  289.     }
  290.    
  291.     treatment<-c(rnormp2(larger.bimodal.proportion*100000,0+effect,treatment.sd1,treatment.shape1),
  292.                  rnormp2(smaller.bimodal.proportion*100000,0+effect+bimodal.offset,treatment.sd2,treatment.shape2))
  293.    
  294.    
  295.       if(batchmode!=T){
  296.         if(using.RStudio!=T){
  297.         dev.new(width=600, height=500)
  298.         }
  299.         }
  300.    
  301.     if(using.RStudio==T){
  302.       live.plot.in.own.window=F
  303.       live.plot=F
  304.     }
  305.     layout(matrix(c(1,2,3,4,5,6), 3, 2, byrow = TRUE))
  306.     hist(treatment, col=rgb(0,0,1,.5), freq=F,
  307.          breaks=seq(min(c(treatment,control))-1,max(c(treatment,control))+1,by=.1),
  308.          xlab="Value", ylim=c(0,1.2),
  309.          main=c("Population Distributions",
  310.                 paste("Difference in Means=", round(effect,2)),
  311.                 paste("Bimodal Proportions=", round(larger.bimodal.proportion,2), ",", round(smaller.bimodal.proportion,2),
  312.                       " Bimodal Offset=", round(bimodal.offset,2))
  313.          ),
  314.          sub=paste("Bimodal Effect on Treatment Only=",bimodal.treatment.only)
  315.     )
  316.     hist(control, col=rgb(1,0,0,.5), freq=F,
  317.          breaks=seq(min(c(treatment,control))-1,max(c(treatment,control))+1,by=.1),
  318.          add=T)
  319.     legend(x="topleft", legend=c(paste("Treatment:",
  320.                                        "mu=", 0+effect,
  321.                                        ", sd1=", round(treatment.sd1,2),
  322.                                        ", p1=", round(treatment.shape1,2),
  323.                                        ", sd2=", round(treatment.sd2,2),
  324.                                        ", p2=", round(treatment.shape2,2)
  325.     ),
  326.                                  paste("Control:",
  327.                                        "mu=",0,
  328.                                        ", sd1=", round(control.sd1,2),
  329.                                        ", p1=", round(control.shape1,2),
  330.                                        ", sd2=", round(control.sd2,2),
  331.                                        ", p2=", round(control.shape2,2))
  332.     ),
  333.            col=c(rgb(0,0,1,.5),rgb(1,0,0,.5)), lwd=6, bty="n")
  334.    
  335.     plot(initial.samp.size,1, type="l", lwd=2, col="White",
  336.          #xlim=c(initial.samp.size,(max.samps+samp.interval)),
  337.          xlim=c(initial.samp.size,(max.samps)),
  338.          xlab="Sample Size (per group)",
  339.          ylim=c(0,1), ylab="P-value",
  340.          main=c(paste("Representative P-value Traces (Significant: P <",pvalue,")"),
  341.                 paste("Outlier Cutoff(",outlier.method,") =",round(outlier.cutoff,2),
  342.                       ", Outlier Bias (High,Low) =", "(", round(outlier.bias.high,2), ",", round(outlier.bias.low,2), ")",
  343.                       sep=""),
  344.                 paste("# of Outcomes=", number.of.outcomes)
  345.          ),
  346.          sub=paste("Sequential Sampling=", sequential, ", U-test=", also.perform.Utest)
  347.     )
  348.     abline(h=pvalue, lwd=2)
  349.    
  350.     ###Start Simulation###
  351.     all.results=NULL
  352.     pb <- txtProgressBar(min = 0, max = runs*max.samps, style = 3)
  353.     for(r in 1:runs){
  354.      
  355.      
  356.      
  357.       ###Start Run###
  358.       run.result=NULL
  359.      
  360.       experiment.results=vector("list", number.of.outcomes)
  361.       for(t in 1:number.of.outcomes){
  362.         if(bimodal.treatment.only==T){
  363.           random.sample<-rep(1,initial.samp.size)
  364.         }else{
  365.           random.sample<-bimodal.pop[runif(initial.samp.size,1,1000)]
  366.         }
  367.         control<-c(rnormp2(length(which(random.sample==1)),0,control.sd1,control.shape1),
  368.                    rnormp2(length(which(random.sample==0)),0+bimodal.offset,control.sd2,control.shape2))
  369.        
  370.         random.sample<-bimodal.pop[runif(initial.samp.size,1,1000)]
  371.         treatment<-c(rnormp2(length(which(random.sample==1)),0+effect,treatment.sd1,treatment.shape1),
  372.                      rnormp2(length(which(random.sample==0)),0+effect+bimodal.offset,treatment.sd2,treatment.shape2))
  373.        
  374.         experiment.result<-rbind(control,treatment)
  375.         experiment.results[[t]]<-experiment.result
  376.       }
  377.      
  378.      
  379.      
  380.       x=initial.samp.size
  381.      
  382.       while(x <(max.samps+samp.interval)){
  383.         sample.result=NULL
  384.         n.cutoff=min(n.cutoff,x)
  385.        
  386.         test.results=NULL
  387.         for(t in 1:number.of.outcomes){
  388.          
  389.           control<-experiment.results[[t]][1,]
  390.           treatment<-experiment.results[[t]][2,]
  391.          
  392.          
  393.           if(also.perform.Utest==T){
  394.             pval1<-t.test(treatment,control)$p.value
  395.             pval2<-wilcox.test(treatment,control)$p.value
  396.             pval<-min(pval1,pval2)
  397.            
  398.             #Drop Outliers
  399.             if(drop.outliers==T){
  400.               if(pval>pvalue){
  401.                 control.temp<-control
  402.                 treatment.temp<-treatment
  403.                
  404.                 if(outlier.method=="SD"){
  405.                 high.cutoff<-mean(control.temp, na.rm=T)+outlier.cutoff*sd(control.temp, na.rm=T)
  406.                 low.cutoff<-mean(control.temp, na.rm=T)-outlier.cutoff*sd(control.temp, na.rm=T)
  407.                 high.outliers<-which(control.temp>high.cutoff)
  408.                 low.outliers<-which(control.temp<low.cutoff)
  409.                
  410.                 control.temp[high.outliers]<-NA
  411.                 control.temp[low.outliers]<-NA
  412.                
  413.                 high.cutoff<-mean(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*sd(treatment.temp, na.rm=T)
  414.                 low.cutoff<-mean(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*sd(treatment.temp, na.rm=T)
  415.                 high.outliers<-which(treatment.temp>high.cutoff)
  416.                 low.outliers<-which(treatment.temp<low.cutoff)
  417.                
  418.                 treatment.temp[high.outliers]<-NA
  419.                 treatment.temp[low.outliers]<-NA
  420.                 }
  421.                
  422.                 if(outlier.method=="MAD"){
  423.                   high.cutoff<-median(control.temp, na.rm=T)+outlier.cutoff*mad(control.temp, na.rm=T)
  424.                   low.cutoff<-median(control.temp, na.rm=T)-outlier.cutoff*mad(control.temp, na.rm=T)
  425.                   high.outliers<-which(control.temp>high.cutoff)
  426.                   low.outliers<-which(control.temp<low.cutoff)
  427.                  
  428.                   control.temp[high.outliers]<-NA
  429.                   control.temp[low.outliers]<-NA
  430.                  
  431.                   high.cutoff<-median(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*mad(treatment.temp, na.rm=T)
  432.                   low.cutoff<-median(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*mad(treatment.temp, na.rm=T)
  433.                   high.outliers<-which(treatment.temp>high.cutoff)
  434.                   low.outliers<-which(treatment.temp<low.cutoff)
  435.                  
  436.                   treatment.temp[high.outliers]<-NA
  437.                   treatment.temp[low.outliers]<-NA
  438.                  
  439.                 }
  440.                            
  441.                
  442.                 if(
  443.                   length(treatment.temp[which(!is.na(treatment.temp))]) > (n.cutoff - 1) &&
  444.                     length(control.temp[which(!is.na(control.temp))]) > (n.cutoff - 1)
  445.                 ){
  446.                  
  447.                   pval1<-t.test(treatment.temp,control.temp)$p.value
  448.                   pval2<-wilcox.test(treatment.temp,control.temp)$p.value
  449.                   pval<-min(pval1,pval2)
  450.                  
  451.                   control<-control.temp
  452.                   treatment<-treatment.temp
  453.                  
  454.                   #Replace Outliers for Small N
  455.                 }else{
  456.                  
  457.                   while(
  458.                     length(treatment.temp[which(!is.na(treatment.temp))]) < n.cutoff &&
  459.                       length(control.temp[which(!is.na(control.temp))]) < n.cutoff
  460.                   ){
  461.                     if(bimodal.treatment.only==T){
  462.                       random.sample<-rep(1,length(which(is.na(control.temp))))
  463.                     }else{
  464.                       random.sample<-bimodal.pop[runif(length(which(is.na(control.temp))),1,1000)]
  465.                     }
  466.                     control.temp[which(is.na(control.temp))]<-c(rnormp2(length(which(random.sample==1)),
  467.                                                                         0,control.sd1,control.shape1),
  468.                                                                 rnormp2(length(which(random.sample==0)),
  469.                                                                         0+bimodal.offset,control.sd2,control.shape2))
  470.                    
  471.                     random.sample<-bimodal.pop[runif(length(which(is.na(treatment.temp))),1,1000)]
  472.                     treatment.temp[which(is.na(treatment.temp))]<-c(rnormp2(length(which(random.sample==1)),
  473.                                                                             0+effect,treatment.sd1,treatment.shape1),
  474.                                                                     rnormp2(length(which(random.sample==0)),
  475.                                                                             0+effect+bimodal.offset,treatment.sd2,treatment.shape2))
  476.                    
  477.                     if(outlier.method=="SD"){
  478.                       high.cutoff<-mean(control.temp, na.rm=T)+outlier.cutoff*sd(control.temp, na.rm=T)
  479.                       low.cutoff<-mean(control.temp, na.rm=T)-outlier.cutoff*sd(control.temp, na.rm=T)
  480.                       high.outliers<-which(control.temp>high.cutoff)
  481.                       low.outliers<-which(control.temp<low.cutoff)
  482.                      
  483.                       control.temp[high.outliers]<-NA
  484.                       control.temp[low.outliers]<-NA
  485.                      
  486.                       high.cutoff<-mean(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*sd(treatment.temp, na.rm=T)
  487.                       low.cutoff<-mean(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*sd(treatment.temp, na.rm=T)
  488.                       high.outliers<-which(treatment.temp>high.cutoff)
  489.                       low.outliers<-which(treatment.temp<low.cutoff)
  490.                      
  491.                       treatment.temp[high.outliers]<-NA
  492.                       treatment.temp[low.outliers]<-NA
  493.                     }
  494.                    
  495.                     if(outlier.method=="MAD"){
  496.                       high.cutoff<-median(control.temp, na.rm=T)+outlier.cutoff*mad(control.temp, na.rm=T)
  497.                       low.cutoff<-median(control.temp, na.rm=T)-outlier.cutoff*mad(control.temp, na.rm=T)
  498.                       high.outliers<-which(control.temp>high.cutoff)
  499.                       low.outliers<-which(control.temp<low.cutoff)
  500.                      
  501.                       control.temp[high.outliers]<-NA
  502.                       control.temp[low.outliers]<-NA
  503.                      
  504.                       high.cutoff<-median(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*mad(treatment.temp, na.rm=T)
  505.                       low.cutoff<-median(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*mad(treatment.temp, na.rm=T)
  506.                       high.outliers<-which(treatment.temp>high.cutoff)
  507.                       low.outliers<-which(treatment.temp<low.cutoff)
  508.                      
  509.                       treatment.temp[high.outliers]<-NA
  510.                       treatment.temp[low.outliers]<-NA
  511.                      
  512.                     }
  513.                    
  514.                   } #end while sample size too small loop
  515.                  
  516.                   pval1<-t.test(treatment.temp,control.temp)$p.value
  517.                   pval2<-wilcox.test(treatment.temp,control.temp)$p.value
  518.                   pval<-min(pval1,pval2)
  519.                   control<-control.temp
  520.                   treatment<-treatment.temp
  521.                  
  522.                 } # end replace pvals
  523.                
  524.               } # end if pval...
  525.             } #end drop outliers
  526.            
  527.             # Do if U-test==F
  528.           }else{
  529.             pval<-t.test(treatment,control)$p.value
  530.            
  531.             #Drop Outliers
  532.             if(drop.outliers==T){
  533.               if(pval>pvalue){
  534.                 control.temp<-control
  535.                 treatment.temp<-treatment
  536.                
  537.                 if(outlier.method=="SD"){
  538.                   high.cutoff<-mean(control.temp, na.rm=T)+outlier.cutoff*sd(control.temp, na.rm=T)
  539.                   low.cutoff<-mean(control.temp, na.rm=T)-outlier.cutoff*sd(control.temp, na.rm=T)
  540.                   high.outliers<-which(control.temp>high.cutoff)
  541.                   low.outliers<-which(control.temp<low.cutoff)
  542.                  
  543.                   control.temp[high.outliers]<-NA
  544.                   control.temp[low.outliers]<-NA
  545.                  
  546.                   high.cutoff<-mean(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*sd(treatment.temp, na.rm=T)
  547.                   low.cutoff<-mean(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*sd(treatment.temp, na.rm=T)
  548.                   high.outliers<-which(treatment.temp>high.cutoff)
  549.                   low.outliers<-which(treatment.temp<low.cutoff)
  550.                  
  551.                   treatment.temp[high.outliers]<-NA
  552.                   treatment.temp[low.outliers]<-NA
  553.                 }
  554.                
  555.                 if(outlier.method=="MAD"){
  556.                   high.cutoff<-median(control.temp, na.rm=T)+outlier.cutoff*mad(control.temp, na.rm=T)
  557.                   low.cutoff<-median(control.temp, na.rm=T)-outlier.cutoff*mad(control.temp, na.rm=T)
  558.                   high.outliers<-which(control.temp>high.cutoff)
  559.                   low.outliers<-which(control.temp<low.cutoff)
  560.                  
  561.                   control.temp[high.outliers]<-NA
  562.                   control.temp[low.outliers]<-NA
  563.                  
  564.                   high.cutoff<-median(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*mad(treatment.temp, na.rm=T)
  565.                   low.cutoff<-median(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*mad(treatment.temp, na.rm=T)
  566.                   high.outliers<-which(treatment.temp>high.cutoff)
  567.                   low.outliers<-which(treatment.temp<low.cutoff)
  568.                  
  569.                   treatment.temp[high.outliers]<-NA
  570.                   treatment.temp[low.outliers]<-NA
  571.                  
  572.                 }
  573.                 if(
  574.                   length(treatment.temp[which(!is.na(treatment.temp))]) > (n.cutoff - 1) &&
  575.                     length(control.temp[which(!is.na(control.temp))]) > (n.cutoff - 1)
  576.                 ){
  577.                   pval<-t.test(control.temp,treatment.temp)$p.value
  578.                   control<-control.temp
  579.                   treatment<-treatment.temp
  580.                
  581.                 #Replace Outliers for Small N
  582.               }else{
  583.                
  584.                 while(
  585.                   length(treatment.temp[which(!is.na(treatment.temp))]) < n.cutoff &&
  586.                     length(control.temp[which(!is.na(control.temp))]) < n.cutoff
  587.                 ){
  588.                   if(bimodal.treatment.only==T){
  589.                     random.sample<-rep(1,length(which(is.na(control.temp))))
  590.                   }else{
  591.                     random.sample<-bimodal.pop[runif(length(which(is.na(control.temp))),1,1000)]
  592.                   }
  593.                   control.temp[which(is.na(control.temp))]<-c(rnormp2(length(which(random.sample==1)),
  594.                                                                       0,control.sd1,control.shape1),
  595.                                                               rnormp2(length(which(random.sample==0)),
  596.                                                                       0+bimodal.offset,control.sd2,control.shape2))
  597.                  
  598.                   random.sample<-bimodal.pop[runif(length(which(is.na(treatment.temp))),1,1000)]
  599.                   treatment.temp[which(is.na(treatment.temp))]<-c(rnormp2(length(which(random.sample==1)),
  600.                                                                           0+effect,treatment.sd1,treatment.shape1),
  601.                                                                   rnormp2(length(which(random.sample==0)),
  602.                                                                           0+effect+bimodal.offset,treatment.sd2,treatment.shape2))
  603.                  
  604.                   if(outlier.method=="SD"){
  605.                     high.cutoff<-mean(control.temp, na.rm=T)+outlier.cutoff*sd(control.temp, na.rm=T)
  606.                     low.cutoff<-mean(control.temp, na.rm=T)-outlier.cutoff*sd(control.temp, na.rm=T)
  607.                     high.outliers<-which(control.temp>high.cutoff)
  608.                     low.outliers<-which(control.temp<low.cutoff)
  609.                    
  610.                     control.temp[high.outliers]<-NA
  611.                     control.temp[low.outliers]<-NA
  612.                    
  613.                     high.cutoff<-mean(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*sd(treatment.temp, na.rm=T)
  614.                     low.cutoff<-mean(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*sd(treatment.temp, na.rm=T)
  615.                     high.outliers<-which(treatment.temp>high.cutoff)
  616.                     low.outliers<-which(treatment.temp<low.cutoff)
  617.                    
  618.                     treatment.temp[high.outliers]<-NA
  619.                     treatment.temp[low.outliers]<-NA
  620.                   }
  621.                  
  622.                   if(outlier.method=="MAD"){
  623.                     high.cutoff<-median(control.temp, na.rm=T)+outlier.cutoff*mad(control.temp, na.rm=T)
  624.                     low.cutoff<-median(control.temp, na.rm=T)-outlier.cutoff*mad(control.temp, na.rm=T)
  625.                     high.outliers<-which(control.temp>high.cutoff)
  626.                     low.outliers<-which(control.temp<low.cutoff)
  627.                    
  628.                     control.temp[high.outliers]<-NA
  629.                     control.temp[low.outliers]<-NA
  630.                    
  631.                     high.cutoff<-median(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*mad(treatment.temp, na.rm=T)
  632.                     low.cutoff<-median(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*mad(treatment.temp, na.rm=T)
  633.                     high.outliers<-which(treatment.temp>high.cutoff)
  634.                     low.outliers<-which(treatment.temp<low.cutoff)
  635.                    
  636.                     treatment.temp[high.outliers]<-NA
  637.                     treatment.temp[low.outliers]<-NA
  638.                    
  639.                   }
  640.                  
  641.                 } #end while sample size too small loop
  642.                
  643.                 pval<-t.test(treatment.temp,control.temp)$p.value
  644.                 control<-control.temp
  645.                 treatment<-treatment.temp
  646.                
  647.               } # end replace pvals
  648.              
  649.              
  650.             }# end if pval>pvalue
  651.           } #end drop outliers
  652.          
  653.          
  654.          
  655.         }#end If U-test==F
  656.        
  657.         measured.effect=mean(treatment, na.rm=T)-mean(control, na.rm=T)
  658.        
  659.         test.result<-cbind(pval,measured.effect)
  660.         test.results<-rbind(test.results,test.result)
  661.       }
  662.      
  663.       pval<-min(test.results[,1])
  664.       measured.effect<-test.results[which(test.results[,1]==min(test.results[,1]))[1],2]
  665.      
  666.      
  667.      
  668.       sample.result<-cbind(sample.result,r, x,pval,measured.effect)
  669.       run.result<-rbind(run.result,sample.result)
  670.      
  671.       ##Update Live Plot
  672.       if(live.plot==T){
  673.         lines(run.result[,2],run.result[,3], type="l", lwd=3, col=rainbow(runs)[r],
  674.               xlim=c(initial.samp.size,max.samps+initial.samp.size),ylim=c(0,1)
  675.         )
  676.        
  677.         if(live.plot.in.own.window==T && x==initial.samp.size && r==1){
  678.           dev.new()
  679.           plot(initial.samp.size,1, type="l", lwd=2, col="White",
  680.                xlim=c(initial.samp.size,(max.samps+samp.interval)), xlab="Sample Size (per group)",
  681.                ylim=c(0,1), ylab="P-value",
  682.                main=c(paste("Representative P-value Traces (Significant: P <",pvalue,")"),
  683.                       paste("Outlier Cutoff(SDs) =",round(outlier.cutoff,2),
  684.                             ", Outlier Bias (High,Low) =", "(", round(outlier.bias.high,2), ",", round(outlier.bias.low,2), ")"),
  685.                       paste("# of Outcomes=", number.of.outcomes)
  686.                ),
  687.                sub=paste("Sequential Sampling=", sequential, ", U-test=", also.perform.Utest)
  688.           )
  689.           abline(h=pvalue, lwd=2)
  690.           lines(run.result[,2],run.result[,3], type="l", lwd=3, col=rainbow(runs)[r],
  691.                 xlim=c(initial.samp.size,max.samps+initial.samp.size),ylim=c(0,1)
  692.           )
  693.           dev.set(which=dev.prev())
  694.         }
  695.        
  696.         if(live.plot.in.own.window==T && x*r>initial.samp.size){
  697.           dev.set(which=dev.next())
  698.           lines(run.result[,2],run.result[,3], type="l", lwd=3, col=rainbow(runs)[r],
  699.                 xlim=c(initial.samp.size,max.samps+initial.samp.size),ylim=c(0,1)
  700.           )
  701.           dev.set(which=dev.prev())
  702.         }
  703.         # End Update Live Plot
  704.        
  705.        
  706.        
  707.         Sys.sleep(sleep.time)
  708.       }
  709.      
  710.      
  711.       #Generate New Samples
  712.       if(sequential==T){
  713.        
  714.         ##Take Sequential Samples##
  715.        
  716.         for(t in 1:number.of.outcomes){
  717.          
  718.           #Control
  719.           if(bimodal.treatment.only==T){
  720.             random.sample<-rep(1,samp.interval)
  721.           }else{
  722.             random.sample<-bimodal.pop[runif(samp.interval,1,1000)]
  723.           }
  724.           control<-append(experiment.results[[t]][1,],c(rnormp2(length(which(random.sample==1)),
  725.                                                                 0,control.sd1,control.shape1),
  726.                                                         rnormp2(length(which(random.sample==0)),
  727.                                                                 0+bimodal.offset,control.sd2,control.shape2)
  728.           ))
  729.          
  730.           #Treatment
  731.           random.sample<-bimodal.pop[runif(samp.interval,1,1000)]
  732.           treatment<-append(experiment.results[[t]][2,],c(rnormp2(length(which(random.sample==1)),
  733.                                                                   0+effect,treatment.sd1,treatment.shape1),
  734.                                                           rnormp2(length(which(random.sample==0)),
  735.                                                                   0+effect+bimodal.offset,treatment.sd2,treatment.shape2)
  736.           ))
  737.          
  738.           experiment.result<-rbind(control,treatment)
  739.           experiment.results[[t]]<-experiment.result
  740.          
  741.         }#End Sequential
  742.        
  743.       }else{
  744.        
  745.         ##Take Independant Samples (non sequential)##
  746.         experiment.results=vector("list", number.of.outcomes)
  747.         for(t in 1:number.of.outcomes){
  748.          
  749.           #Control
  750.           if(bimodal.treatment.only==T){
  751.             random.sample<-rep(1,x)
  752.           }else{
  753.             random.sample<-bimodal.pop[runif(x,1,1000)]
  754.           }
  755.           control<-c(rnormp2(length(which(random.sample==1)),
  756.                              0,control.sd1,control.shape1),
  757.                      rnormp2(length(which(random.sample==0)),
  758.                              0+bimodal.offset,control.sd2,control.shape2)
  759.           )
  760.          
  761.           #Treatment
  762.           random.sample<-bimodal.pop[runif(x,1,1000)]
  763.           treatment<-c(rnormp2(length(which(random.sample==1)),
  764.                                0+effect,treatment.sd1,treatment.shape1),
  765.                        rnormp2(length(which(random.sample==0)),
  766.                                0+effect+bimodal.offset,treatment.sd2,treatment.shape2)
  767.           )
  768.          
  769.           experiment.result<-rbind(control,treatment)
  770.           experiment.results[[t]]<-experiment.result
  771.          
  772.         }       #End independant
  773.        
  774.       } #End sampling
  775.      
  776.      
  777.       setTxtProgressBar(pb, x + max.samps*(r-1))
  778.       x<-x+samp.interval
  779.      
  780.     }   # End Run
  781.     all.results<-rbind(all.results,run.result)
  782.    
  783.     if(live.plot==F && r==runs){
  784.       colnum<-1
  785.       for(k in (runs-9):runs){
  786.         xvar<-all.results[which(all.results[,1]==k),2]
  787.         yvar<-all.results[which(all.results[,1]==k),3]
  788.         lines(xvar,yvar, type="l", lwd=3, col=rainbow(10)[colnum],
  789.               xlim=c(initial.samp.size,max.samps+initial.samp.size),ylim=c(0,1)
  790.         )
  791.         colnum<-colnum+1
  792.       }
  793.     }
  794.    
  795.   }     # End Simulation
  796.   close(pb)
  797.  
  798.  
  799.  
  800.   ###Extract Significant Results###
  801.   rownames(all.results)<-seq(1,nrow(all.results))
  802.   sigs<-which(all.results[,3]<pvalue)
  803.   percent.sigs<-100*hist(all.results[sigs,2], plot=F,
  804.                          breaks=seq(initial.samp.size-1, max.samps+2*samp.interval,by=(samp.interval))
  805.   )$counts/runs
  806.  
  807.   ###Extract First Significant Result for each Simulation Run###
  808.   first.sigs=NULL
  809.   first.sigs.filtered=NULL
  810.   for(i in 1:runs){
  811.     temp<-all.results[which(all.results[,1]==i),]
  812.     if(length(temp[which(temp[,3]<pvalue),2])==0){
  813.       first.sig<-Inf
  814.     }else{
  815.       first.sig<-min(temp[which(temp[,3]<pvalue),2])
  816.     }
  817.     first.sigs<-rbind(first.sigs,first.sig)
  818.   }
  819.  
  820.   first.sigs.filtered<-first.sigs[which(first.sigs!=Inf),]
  821.  
  822.   percent.atLeastOne.sig<-100*cumsum(
  823.     hist(first.sigs.filtered,
  824.          breaks=seq(initial.samp.size-1, max.samps+2*samp.interval,by=(samp.interval)),
  825.          plot=F
  826.     )$counts
  827.   )/runs
  828.  
  829.  
  830.   ###Create Summary Plots###
  831.  
  832.  
  833.   if(bimodal.treatment.only==T){
  834.    
  835.     plot(all.results[sigs,2],all.results[sigs,4],
  836.          xlab="Sample Size (per group)",xlim=c(0, max.samps+samp.interval),
  837.          ylab="Measured Effect",
  838.          ylim=c(min(c(-1,all.results[sigs,4]),na.rm=T),max(c(1,all.results[sigs,4]),na.rm=T)),
  839.          main=c("Measured Effect Sizes of Significant Results",
  840.                 paste("Significance Level: P <", pvalue),
  841.                 paste("True Difference=", round(effect,2), "+",
  842.                       round(bimodal.offset,2), "for", round(100*smaller.bimodal.proportion,2), "% of Treatment Group")),
  843.          col="Cyan")
  844.     abline(h=effect, lwd=2)
  845.     abline(h=bimodal.offset, lwd=1, lty=2)
  846.    
  847.    
  848.   }else{
  849.    
  850.     plot(all.results[sigs,2],all.results[sigs,4],
  851.          xlab="Sample Size (per group)",xlim=c(0, max.samps+samp.interval),
  852.          ylab="Measured Effect", ylim=c(min(c(-1,all.results[sigs,4]),na.rm=T),max(c(1,all.results[sigs,4]),na.rm=T)),
  853.          main=c("Measured Effect Sizes of Significant Results",
  854.                 paste("Significance Level: P <", pvalue),
  855.                 paste("True Difference=", effect)),
  856.          col="Cyan")
  857.     abline(h=effect, lwd=2)
  858.   }
  859.  
  860.     pval.max.percents=NULL
  861.     for(s in 1:length(unique(all.results[,2]))){
  862. size<-unique(all.results[,2])[s]
  863.       temp<-max(100*hist(all.results[which(all.results[,2]==size),3],
  864.                          breaks=seq(0,1,by=.005),
  865.                          plot=F)$counts/runs)
  866.       pval.max.percents<-rbind(pval.max.percents,temp)
  867.     }
  868.     ymax<-min(1.5*max(pval.max.percents),100)
  869.    
  870.     Sys.sleep(.1)
  871.     plot(seq(0,.995,by=.005),
  872.          100*hist(all.results[which(all.results[,2]==initial.samp.size),3],
  873.                   breaks=seq(0,1,by=.005),
  874.                   plot=F)$counts/runs,
  875.          ylab="Percent of P-values", ylim=c(0,ymax),
  876.          xlab="P-value", xlim=c(0,1),
  877.          type="l", col=rainbow(length(unique(all.results[,2])))[1],
  878.          main="Distribution of All P-values"
  879.     )
  880.     for(s in 2:length(unique(all.results[,2]))){
  881.       size<-unique(all.results[,2])[s]
  882.       lines(seq(0,.995,by=.005),
  883.             100*hist(all.results[which(all.results[,2]==size),3],
  884.                      breaks=seq(0,1,by=.005),
  885.                      plot=F)$counts/runs,
  886.             type="l",col=rainbow(length(unique(all.results[,2])))[s]
  887.       )
  888.      
  889.     }
  890.     abline(h=100/length(seq(0,.995,by=.005)))
  891.    
  892.     image.plot(matrix(seq(initial.samp.size,max.samps,by=samp.interval)),
  893.                col=rainbow(length(unique(all.results[,2]))),
  894.                smallplot=c(.6,.85,.7,.75),
  895.                legend.only=T, horizontal=T, legend.shrink=.5, add=T,
  896.                legend.args=list(text="Sample Size", side=3, line=.2, font=2, cex=.5, bty="o", bg="White"),
  897.                axis.args=list(
  898.                  at=c(initial.samp.size,max.samps),
  899.                               labels=as.character(c(initial.samp.size,max.samps
  900.                               )
  901.                               )
  902.                )
  903.     )
  904.  
  905.   xsamps<-unique(all.results[,2])
  906.  
  907.   if(sequential==T){
  908.     plot(xsamps,
  909.          percent.atLeastOne.sig[1:length(xsamps)],
  910.          ylab=expression("Pr(# Significant)">=1), ylim=c(0,max(percent.atLeastOne.sig)+10),
  911.          xlab="Sample Size (per Group)",xlim=c(0, max.samps+samp.interval),
  912.          type="s", col="Cyan", lwd=5,
  913.          main=c("Probability of at Least 1 Significant Result",
  914.                 paste("Number of Simulations=",runs),
  915.                 paste("Initial Sample Size=",initial.samp.size, " Sampling Interval=",samp.interval)
  916.          )
  917.     )
  918.     abline(h=seq(0,max(percent.atLeastOne.sig)+10,by=10),
  919.            lwd=1, col=rgb(0,0,0,.5))
  920.     abline(v=seq(0,max.samps+initial.samp.size,by=10), lwd=1, col=rgb(0,0,0,.5))
  921.     abline(h=pvalue*100,
  922.            lwd=2, col=rgb(0,0,0,1))
  923.    
  924.   }else{
  925.     plot(xsamps,
  926.          percent.sigs[1:length(xsamps)],
  927.          ylim=c(0,max(percent.sigs)+10), ylab="% Significant Results",
  928.          xlim=c(0, max.samps+samp.interval), xlab="Sample Size (per Group)",
  929.          type="l", lwd=4, col="Cyan",
  930.          main=c("Percent of Significant Results for Different Sample Sizes",
  931.                 paste("Number of Simulations=",runs),
  932.                 paste("Initial Sample Size=",initial.samp.size, " Sampling Interval=",samp.interval)
  933.          )
  934.     )
  935.     abline(h=seq(0,max(percent.sigs)+10,by=10),
  936.            lwd=1, col=rgb(0,0,0,.5))
  937.     abline(v=seq(0,max.samps+initial.samp.size,by=10), lwd=1, col=rgb(0,0,0,.5))
  938.     abline(h=pvalue*100, lwd=2)
  939.   }
  940.     Sys.sleep(.1)
  941.     plot(seq(0,.995,by=.005),
  942.          100*hist(all.results[which(all.results[,2]==initial.samp.size),3],
  943.                   breaks=seq(0,1,by=.005),
  944.                   plot=F)$counts/runs,
  945.          ylab="Percent of P-values", ylim=c(0,ymax),
  946.          xlab="P-value", xlim=c(0,pvalue),
  947.          type="l", col=rainbow(length(unique(all.results[,2])))[1],
  948.          main="Distribution of Significant P-values"
  949.     )
  950.     for(s in 2:length(unique(all.results[,2]))){
  951.       size<-unique(all.results[,2])[s]
  952.       lines(seq(0,.995,by=.005),
  953.             100*hist(all.results[which(all.results[,2]==size),3],
  954.                      breaks=seq(0,1,by=.005),
  955.                      plot=F)$counts/runs,
  956.             type="l",col=rainbow(length(unique(all.results[,2])))[s]
  957.       )
  958.      
  959.     }
  960.     abline(h=100/length(seq(0,.995,by=.005)))
  961.    
  962.     image.plot(matrix(seq(initial.samp.size,max.samps,by=samp.interval)),
  963.                col=rainbow(length(unique(all.results[,2]))),
  964.                smallplot=c(.6,.85,.7,.75),
  965.                legend.only=T, horizontal=T, legend.shrink=.3, add=T,
  966.                legend.args=list(text="Sample Size", side=3, line=.2, font=2, cex=.5),
  967.                axis.args=list(
  968.                  at=c(initial.samp.size,max.samps),
  969.                  labels=as.character(c(initial.samp.size,max.samps
  970.                  )
  971.                  )
  972.                )
  973.     )
  974.    
  975.    
  976. } #End Script
  977.  
  978. if(batchmode==T){
  979.   effect.info=NULL
  980.   for(samp.point in 1:length(sort(unique(all.results[sigs,2])))){
  981.     sp<-sort(unique(all.results[sigs,2]))[samp.point]
  982.     temp<-all.results[sigs,4][which(all.results[sigs,2]==sp)]
  983.     effect.temp2<-cbind(
  984.       max.pos<-max(temp[which(temp>0)]),
  985.       mean.pos<-mean(temp[which(temp>0)]),
  986.       min.pos<-min(temp[which(temp>0)]),
  987.       max.neg<-max(temp[which(temp<0)]),
  988.       mean.neg<-mean(temp[which(temp<0)]),
  989.       min.neg<-min(temp[which(temp<0)])
  990.     )
  991.    
  992.     effect.info<-rbind(effect.info,effect.temp2)
  993.   }# End Batch
  994.  
  995.   batch.result<-cbind(b,
  996.                       sort(unique(all.results[sigs,2])),
  997.                       percent.sigs[1:length(xsamps)],
  998.                       percent.atLeastOne.sig[1:length(xsamps)],
  999.                       effect.info
  1000.   )
  1001.   colnames(batch.result)<-c("Batch", "Samp.Size", "Percent.Sig",
  1002.                             "Percent.atLeastOne.sig","max.pos","mean.pos","min.pos",
  1003.                             "max.neg","mean.neg","min.neg")
  1004.  
  1005.   batch.results<-rbind(batch.results,batch.result)
  1006.  
  1007.   sim.params<-cbind(rep(b,12),
  1008.                     rbind(
  1009.                       treatment.sd1,
  1010.                       treatment.shape1,
  1011.                       treatment.sd2,
  1012.                       treatment.shape2,
  1013.                       control.sd1,
  1014.                       control.shape1,
  1015.                       control.sd2,
  1016.                       control.shape2,
  1017.                       larger.bimodal.proportion,
  1018.                       bimodal.offset,
  1019.                       bimodal.treatment.only,
  1020.                       sequential,      
  1021.                       also.perform.Utest,
  1022.                       number.of.outcomes
  1023.                     )
  1024.   )
  1025.   batch.params<-rbind(batch.params,sim.params)
  1026. }
  1027.  
  1028. print(paste("Batch #", b, " of", batch.size))
  1029. }
  1030.  
  1031. if(batchmode==T){
  1032.   one.sig=NULL
  1033.   for (i in 1:length(unique(batch.results[,2]))){
  1034.     temp1<-unique(batch.results[,2])[i]
  1035.     temp2<-batch.results[which(batch.results[,2]==temp1),]
  1036.     temp3<-cbind(temp1,
  1037.                  temp2[which(temp2[,4]==max(temp2[,4])),1],
  1038.                  temp2[which(temp2[,4]==max(temp2[,4])),4]
  1039.     )
  1040.     one.sig<-rbind(one.sig,temp3)
  1041.   }
  1042.   colnames(one.sig)<-c("Samp.Size", "Batch", "Percent.atleastOne.Significant")
  1043.   rownames(one.sig)<-1:nrow(one.sig)
  1044.  
  1045.  
  1046.   sampsize.sig=NULL
  1047.   for (i in 1:length(unique(batch.results[,2]))){
  1048.     temp1<-unique(batch.results[,2])[i]
  1049.     temp2<-batch.results[which(batch.results[,2]==temp1),]
  1050.     temp3<-cbind(temp1,
  1051.                  temp2[which(temp2[,3]==max(temp2[,3])),1],
  1052.                  temp2[which(temp2[,3]==max(temp2[,3])),3]
  1053.     )
  1054.     sampsize.sig<-rbind(sampsize.sig,temp3)
  1055.   }
  1056.   colnames(sampsize.sig)<-c("Samp.Size", "Batch", "Percent.Sig")
  1057.   rownames(sampsize.sig)<-1:nrow(sampsize.sig)
  1058.  
  1059.   one.sig
  1060.   sampsize.sig
  1061.  
  1062.   batch.params[which(batch.params[,1]==90),]
  1063.  
  1064. }       # End Batches
  1065.  
  1066.  
  1067.  
  1068.  
  1069.  
  1070. verify=F
  1071. if(verify==T){
  1072.  
  1073.  
  1074.   #control<-c(rnormp2(larger.bimodal.proportion*10000,0,control.sd1,control.shape1),
  1075.   #rnormp2(smaller.bimodal.proportion*10000,0+bimodal.offset,control.sd2,control.shape2))
  1076.  
  1077.  
  1078.   #treatment<-c(rnormp2(larger.bimodal.proportion*10000,0+effect,treatment.sd1,treatment.shape1),
  1079.   #rnormp2(smaller.bimodal.proportion*10000,0+effect+bimodal.offset,treatment.sd2,treatment.shape2))
  1080.  
  1081.  
  1082.   #hist(control)
  1083.  
  1084.  
  1085.   require(sciplot)
  1086.  
  1087.   effect=0.1
  1088.   nsamp=1000
  1089.   drop.outliers=T
  1090.   outlier.cutoff<-2
  1091.   out=NULL
  1092.   for(i in 1:nsamp){
  1093.    
  1094.     dev.new()
  1095.     par(mfrow=c(1,2))
  1096.     #a<-c(rnorm(8,10,1),rnorm(2,12,1))
  1097.     #b<-c(rnorm(8,10,1),rnorm(2,12,1))
  1098.    
  1099.     random.sample<-bimodal.pop[runif(initial.samp.size,1,1000)]
  1100.     random.sample
  1101.     control<-c(rnormp2(length(which(random.sample==1)),0,control.sd1,control.shape1),
  1102.                rnormp2(length(which(random.sample==0)),0+bimodal.offset,control.sd2,control.shape2))
  1103.    
  1104.     random.sample<-bimodal.pop[runif(initial.samp.size,1,1000)]
  1105.     random.sample
  1106.     treatment<-c(rnormp2(length(which(random.sample==1)),0+effect,treatment.sd1,treatment.shape1),
  1107.                  rnormp2(length(which(random.sample==0)),0+effect+bimodal.offset,treatment.sd2,treatment.shape2))
  1108.    
  1109.     pval<-t.test(control,treatment)$p.value
  1110.    
  1111.     if(drop.outliers==T){
  1112.       if(pval>pvalue){
  1113.         control.temp<-control
  1114.         treatment.temp<-treatment
  1115.        
  1116.         high.cutoff<-mean(control.temp)+outlier.cutoff*sd(control.temp)
  1117.         low.cutoff<-mean(control.temp)-outlier.cutoff*sd(control.temp)
  1118.         high.outliers<-which(control.temp>high.cutoff)
  1119.         low.outliers<-which(control.temp<low.cutoff)
  1120.        
  1121.         control.temp[high.outliers]<-NA
  1122.         control.temp[low.outliers]<-NA
  1123.        
  1124.         high.cutoff<-mean(treatment.temp)+outlier.cutoff*sd(treatment.temp)
  1125.         low.cutoff<-mean(treatment.temp)-outlier.cutoff*sd(treatment.temp)
  1126.         high.outliers<-which(treatment.temp>high.cutoff)
  1127.         low.outliers<-which(treatment.temp<low.cutoff)
  1128.        
  1129.         treatment.temp[high.outliers]<-NA
  1130.         treatment.temp[low.outliers]<-NA
  1131.        
  1132.         if(
  1133.           length(treatment.temp[which(!is.na(treatment.temp))])>1 &&
  1134.             length(control.temp[which(!is.na(control.temp))])>1
  1135.         ){
  1136.           pval<-t.test(control.temp,treatment.temp)$p.value
  1137.           control<-control.temp
  1138.           treatment<-treatment.temp
  1139.         }
  1140.        
  1141.       }
  1142.     } #end drop outliers
  1143.    
  1144.     adjust=6
  1145.     a<-control+adjust
  1146.     b<-treatment+adjust
  1147.    
  1148.    
  1149.     boxplot(a,b,col=c("Blue","Red"),ylim=c(0,3*adjust))
  1150.     stripchart(a, add=T, at=1,
  1151.                ylim=c(0,3*adjust),
  1152.                vertical=T, pch=16)
  1153.     stripchart(b, add=T, at=2,
  1154.                vertical=T, pch=16)
  1155.    
  1156.     x.fact<-c(rep("a",length(a)),rep("b",length(b)))
  1157.     vals<-c(a,b)
  1158.     bargraph.CI(x.fact,vals, col=c("Blue","Red"), ylim=c(0,3*adjust),
  1159.                 main=round(t.test(a,b)$p.value,3)
  1160.     )
  1161.     control
  1162.     treatment
  1163.    
  1164.    
  1165.     out<-rbind(out,t.test(a,b)$p.value)
  1166.     Sys.sleep(.5)
  1167.     dev.off()
  1168.   }
  1169.  
  1170.   length(which(out<pvalue))/nsamp
  1171.  
  1172.  
  1173.  
  1174. }
  1175.  
  1176. #one.sig
  1177. #sampsize.sig
  1178.  
  1179. #batch.results
  1180. #batch.params[which(batch.params[,1]==40),]
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top