Always keep the mailing list in cc. The code runs for each row in the data. However I get the feeling that there is a mismatch between what you think that is in the data and the actual data. ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance Kliniekstraat 25 1070 Anderlecht Belgium
To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey 2016-04-18 10:35 GMT+02:00 tan sj <sj_style_1...@outlook.com>: > Thanks but it seem like the problem of looping through data is still the > same....i am really wondering where is the mistake.... > > ________________________________________ > From: Thierry Onkelinx <thierry.onkel...@inbo.be> > Sent: Monday, April 18, 2016 7:21 AM > To: tan sj > Cc: r-help > Subject: Re: [R] R [coding : do not run for every row ] > > You can make this much more readable with apply functions. > > result <- apply( > all_combine1, > 1, > function(x){ > p.value <- sapply( > seq_len(nSims), > function(sim){ > gamma1 <- rgamma(x["m"], x["sp(skewness1.5)"], x["scp1"]) > gamma2 <- rgamma(x["n"], x["scp1"], 1) > gamma1 <- gamma1 - x["sp(skewness1.5)"] * x["scp1"] > gamma2 <- gamma2 - x["sp(skewness1.5)"] > c( > equal = t.test(gamma1, gamma2, var.equal=TRUE)$p.value, > unequal = t.test(gamma1,gamma2,var.equal=FALSE)$p.value, > mann = wilcox.test(gamma1,gamma2)$p.value > ) > } > ) > rowMeans(p.value <= alpha) > } > ) > cbind(all_combine1, t(result)) > ir. Thierry Onkelinx > Instituut voor natuur- en bosonderzoek / Research Institute for Nature > and Forest > team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance > Kliniekstraat 25 > 1070 Anderlecht > Belgium > > To call in the statistician after the experiment is done may be no > more than asking him to perform a post-mortem examination: he may be > able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher > The plural of anecdote is not data. ~ Roger Brinner > The combination of some data and an aching desire for an answer does > not ensure that a reasonable answer can be extracted from a given body > of data. ~ John Tukey > > > 2016-04-18 9:05 GMT+02:00 tan sj <sj_style_1...@outlook.com>: >> Hi, i am sorry, the output should be values between 0 and 0.1 and not >> supposed to be 1.00, it is because they are type 1 error rate. And now i get >> output 1.00 for several samples,rhis is no correct. The loop do not run for >> every row. i do not know where is my mistake. As i use the same concept on >> normal distribution setup, i get the result. >> >> Sent from my phone >> >> On Thierry Onkelinx <thierry.onkel...@inbo.be>, Apr 18, 2016 2:55 PM wrote: >> Dear anonymous, >> >> The big mistake in the output might be obvious to you but not to >> others. Please make clear what the correct output should be or at >> least what is wrong with the current output. >> >> And please DO read the posting guide which asks you not to post in HTML. >> ir. Thierry Onkelinx >> Instituut voor natuur- en bosonderzoek / Research Institute for Nature >> and Forest >> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance >> Kliniekstraat 25 >> 1070 Anderlecht >> Belgium >> >> To call in the statistician after the experiment is done may be no >> more than asking him to perform a post-mortem examination: he may be >> able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher >> The plural of anecdote is not data. ~ Roger Brinner >> The combination of some data and an aching desire for an answer does >> not ensure that a reasonable answer can be extracted from a given body >> of data. ~ John Tukey >> >> >> 2016-04-17 19:59 GMT+02:00 tan sj <sj_style_1...@outlook.com>: >>> i have combined all the variables in a matrix, and i wish to conduct a >>> simulation row by row. >>> >>> But i found out the code only works for the every first row after a cycle >>> of nine samples. >>> >>> But after check out the code, i don know where is my mistake... >>> >>> can anyone pls help .... >>> >>> >>> #For gamma disribution with equal skewness 1.5 >>> >>> #to evaluate the same R function on many different sets of data >>> library(parallel) >>> >>> nSims<-100 >>> alpha<-0.05 >>> >>> #set nrow =nsims because wan storing every p-value simulated >>> #for gamma distribution with equal skewness >>> matrix2_equal <-matrix(0,nrow=nSims,ncol=3) >>> matrix5_unequal<-matrix(0,nrow=nSims,ncol=3) >>> matrix8_mann <-matrix(0,nrow=nSims,ncol=3) >>> >>> # to ensure the reproducity of the result >>> #here we declare the random seed generator >>> set.seed(1) >>> >>> ## Put the samples sizes into matrix then use a loop for sample sizes >>> >>> sample_sizes<-matrix(c(10,10,10,25,25,25,25,50,25,100,50,25,50,100,100,25,100,100), >>> nrow=2) >>> >>> #shape parameter for both gamma distribution for equal skewness >>> #forty five cases for each skewness!! >>> shp<-rep(16/9,each=5) >>> >>> #scale parameter for sample 1 >>> #scale paramter for sample 2 set as constant 1 >>> scp1<-c(1,1.5,2,2.5,3) >>> >>> #get all combinations with one row of the sample_sizes matrix >>> ##(use expand.grid)to create a data frame from combination of data >>> >>> ss_sd1<- expand.grid(sample_sizes[2,],shp) >>> scp1<-rep(scp1,9) >>> >>> std2<-rep(sd2,9) >>> >>> #create a matrix combining the forty five cases of combination of sample >>> sizes,shape and scale parameter >>> all_combine1 <- cbind(rep(sample_sizes[1,], 5),ss_sd1,scp1) >>> >>> # name the column samples 1 and 2 and standard deviation >>> colnames(all_combine1) <- c("m", "n","sp(skewness1.5)","scp1") >>> >>> ##for the samples sizes into matrix then use a loop for sample sizes >>> # this loop steps through the all_combine matrix >>> for(ss in 1:nrow(all_combine1)) >>> { >>> #generate samples from the first column and second column >>> m<-all_combine1[ss,1] >>> n<-all_combine1[ss,2] >>> >>> for (sim in 1:nSims) >>> { >>> #generate 2 random samples from gamma distribution with equal >>> skewness >>> gamma1<-rgamma(m,all_combine1[ss,3],all_combine1[ss,4]) >>> gamma2<-rgamma(n,all_combine1[ss,4],1) >>> >>> # minus the population mean to ensure that there is no lose of >>> equality of mean >>> gamma1<-gamma1-all_combine1[ss,3]*all_combine1[ss,4] >>> gamma2<-gamma2-all_combine1[ss,3] >>> >>> #extract p-value out and store every p-value into matrix >>> matrix2_equal[sim,1]<-t.test(gamma1,gamma2,var.equal=TRUE)$p.value >>> >>> matrix5_unequal[sim,2]<-t.test(gamma1,gamma2,var.equal=FALSE)$p.value >>> matrix8_mann[sim,3] <-wilcox.test(gamma1,gamma2)$p.value >>> } >>> ##store the result >>> equal[ss]<- mean(matrix2_equal[,1]<=alpha) >>> unequal[ss]<-mean(matrix5_unequal[,2]<=alpha) >>> mann[ss]<- mean(matrix8_mann[,3]<=alpha) >>> } >>> >>> g_equal<-cbind(all_combine1, equal, unequal, mann) >>> >>> It is my result but it show a very big mistake ....TT >>> m n sp(skewness1.5) scp1 equal unequal mann >>> 1 10 10 1.777778 1.0 0.36 0.34 0.34 >>> 2 10 25 1.777778 1.5 0.84 0.87 0.90 >>> 3 25 25 1.777778 2.0 1.00 1.00 1.00 >>> 4 25 50 1.777778 2.5 1.00 1.00 1.00 >>> 5 25 100 1.777778 3.0 1.00 1.00 1.00 >>> 6 50 25 1.777778 1.0 0.77 0.77 0.84 >>> 7 50 100 1.777778 1.5 1.00 1.00 1.00 >>> 8 100 25 1.777778 2.0 1.00 1.00 1.00 >>> 9 100 100 1.777778 2.5 1.00 1.00 1.00 >>> 10 10 10 1.777778 3.0 1.00 1.00 1.00 >>> 11 10 25 1.777778 1.0 0.48 0.30 0.55 >>> 12 25 25 1.777778 1.5 0.99 0.99 1.00 >>> 13 25 50 1.777778 2.0 1.00 1.00 1.00 >>> 14 25 100 1.777778 2.5 1.00 1.00 1.00 >>> 15 50 25 1.777778 3.0 1.00 1.00 1.00 >>> 16 50 100 1.777778 1.0 0.97 0.97 1.00 >>> 17 100 25 1.777778 1.5 1.00 1.00 1.00 >>> 18 100 100 1.777778 2.0 1.00 1.00 1.00 >>> 19 10 10 1.777778 2.5 1.00 1.00 1.00 >>> 20 10 25 1.777778 3.0 1.00 1.00 1.00 >>> 21 25 25 1.777778 1.0 0.63 0.63 0.71 >>> 22 25 50 1.777778 1.5 0.99 0.99 0.99 >>> 23 25 100 1.777778 2.0 1.00 1.00 1.00 >>> 24 50 25 1.777778 2.5 1.00 1.00 1.00 >>> 25 50 100 1.777778 3.0 1.00 1.00 1.00 >>> 26 100 25 1.777778 1.0 0.83 0.90 0.88 >>> 27 100 100 1.777778 1.5 1.00 1.00 1.00 >>> 28 10 10 1.777778 2.0 1.00 1.00 1.00 >>> 29 10 25 1.777778 2.5 1.00 1.00 1.00 >>> 30 25 25 1.777778 3.0 1.00 1.00 1.00 >>> 31 25 50 1.777778 1.0 0.71 0.66 0.81 >>> 32 25 100 1.777778 1.5 1.00 1.00 1.00 >>> 33 50 25 1.777778 2.0 1.00 1.00 1.00 >>> 34 50 100 1.777778 2.5 1.00 1.00 1.00 >>> 35 100 25 1.777778 3.0 1.00 1.00 1.00 >>> 36 100 100 1.777778 1.0 0.99 0.99 1.00 >>> 37 10 10 1.777778 1.5 0.65 0.65 0.71 >>> 38 10 25 1.777778 2.0 1.00 1.00 1.00 >>> 39 25 25 1.777778 2.5 1.00 1.00 1.00 >>> 40 25 50 1.777778 3.0 1.00 1.00 1.00 >>> 41 25 100 1.777778 1.0 0.90 0.89 0.96 >>> 42 50 25 1.777778 1.5 0.99 0.99 1.00 >>> 43 50 100 1.777778 2.0 1.00 1.00 1.00 >>> 44 100 25 1.777778 2.5 1.00 1.00 1.00 >>> 45 100 100 1.777778 3.0 1.00 1.00 1.00 >>>> >>> >>> >>> >>> [[alternative HTML version deleted]] >>> >>> ______________________________________________ >>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see >>> https://stat.ethz.ch/mailman/listinfo/r-help >>> PLEASE do read the posting guide >>> http://www.R-project.org/posting-guide.html >>> and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.