[R] R [coding : do not run for every row ]

Thierry Onkelinx thierry.onkelinx at inbo.be
Mon Apr 18 12:03:04 CEST 2016


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_1125 op 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.onkelinx op 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_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