API tools faq
paste
Advertisement
Guest User

R simulation

a guest
Jan 30th, 2013
259
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 45.95 KB | None | 0 0
raw download clone embed print report
  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),]
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!