[R] R [coding : do not run for every row ]
Thierry Onkelinx
thierry.onkelinx at inbo.be
Mon Apr 18 09:21:43 CEST 2016
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_1125 op 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.onkelinx op 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_1125 op 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 op 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.
More information about the R-help
mailing list