Permutation Test Example PDF

> #permutation_test_example.R > #perform a permutation test for H0:mu1=mu2 versus H1:mu1 NE mu2 > #from two independent samples from populations with equal variance > > library(psych) #includes the describe function > > setwd("C:/Users/jmard/OneDrive/Data101_Fall2024/Output") > #save gra...

> #permutation_test_example.R > #perform a permutation test for H0:mu1=mu2 versus H1:mu1 NE mu2 > #from two independent samples from populations with equal variance > > library(psych) #includes the describe function > > setwd("C:/Users/jmard/OneDrive/Data101_Fall2024/Output") > #save graphics output as a pdf - saves graph(s) in working directory > pdf(file="permutation_test_example_out.pdf") > > #What are p-values? > #p-values are the probability of obtaining a test statistic at least > #as extreme as the test statistic computed from sample data, > #assuming the null hypothesis is true (mu1=mu2). > > #Permutation tests sample from an estimated null distribution for the data, > #and use the samples to estimate p-values for hypothesis tests. > # > #Simple Example: > #Data: Two independent samples of rocket propellant type (typeA, typeB) > #The response is burn time > > typeA = c(65,81,57,66,82,82,67,59,75,70) #n1=10 > typeB = c(64,71,83,59,65,56,69,74,82,79) #n2=10 > > #generate descriptive statistics > > describe(typeA) vars n mean sd median trimmed mad min max range skew kurtosis se X1 1 10 70.4 9.26 68.5 70.62 11.86 57 82 25 0.03 -1.65 2.93 > describe(typeB) vars n mean sd median trimmed mad min max range skew kurtosis se X1 1 10 70.2 9.37 70 70.38 11.12 56 83 27 -0.02 -1.56 2.96 > > #perform the two-sample t-test > #which is a competitor to the permutation test > > t.test(typeA,typeB,var.equal=TRUE, paired=FALSE) Two Sample t-test data: typeA and typeB t = 0.048008, df = 18, p-value = 0.9622 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -8.552441 8.952441 sample estimates: mean of x mean of y 70.4 70.2 > > #critical value alpha=0.05 > qt(0.975, 18) # 97.5% quantile from t-distribution with 18 df 2.100922 > > #this next line computes the two-tail p-value > pt(0.048008,18,lower.tail = FALSE) + pt(-0.048008,18,lower.tail = TRUE) 0.9622385 > > #Permutation tests sample from an estimated null distribution for the data, > #and use the samples to estimate p-values for hypothesis tests. > # > #Simple Example: > #Data: Two independent samples of rocket propellant type (typeA, typeB) > #The response is burn time > > typeA = c(65,81,57,66,82,82,67,59,75,70) #n1=10 > typeB = c(64,71,83,59,65,56,69,74,82,79) #n2=10 > > #generate descriptive statistics > > describe(typeA) vars n mean sd median trimmed mad min max range skew kurtosis se X1 1 10 70.4 9.26 68.5 70.62 11.86 57 82 25 0.03 -1.65 2.93 > describe(typeB) vars n mean sd median trimmed mad min max range skew kurtosis se X1 1 10 70.2 9.37 70 70.38 11.12 56 83 27 -0.02 -1.56 2.96 > > #perform the two-sample t-test > #which is a competitor to the permutation test > > t.test(typeA,typeB,var.equal=TRUE, paired=FALSE) Two Sample t-test data: typeA and typeB t = 0.048008, df = 18, p-value = 0.9622 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -8.552441 8.952441 sample estimates: mean of x mean of y 70.4 70.2 > > qqnorm(typeA,datax=TRUE) #graphical test of normality > qqline(typeA,datax=TRUE) > > qqnorm(typeB,datax=TRUE) #graphical test of normality > qqline(typeB,datax=TRUE) > > > #generate a single column vector of two variables, type and burntime > example #print data as a column > example type burntime 1 A 65 2 A 81 3 A 57 4 A 66 5 A 82 6 A 82 7 A 67 8 A 59 9 A 75 10 A 70 11 B 64 12 B 71 13 B 83 14 B 59 15 B 65 16 B 56 17 B 69 18 B 74 19 B 82 20 B 79 > > t.test(example$burntime~example$type,var.equal=TRUE) #assuming equal variance Two Sample t-test data: example$burntime by example$type t = 0.048008, df = 18, p-value = 0.9622 alternative hypothesis: true difference in means between group A and group B is not equal to 0 95 percent confidence interval: -8.552441 8.952441 sample estimates: mean in group A mean in group B 70.4 70.2 > > if (FALSE) + {" + The logic behind a permutation test is straightforward. These are the 20 (in this example) scores we observed. + If in fact the two samples came from the same population (assume the null hypothesis) + then here is one possible permutation of the twenty scores. There are, in fact 20 Choose 10 permutations + "} > > choose(20,10) # 184756 =20!/(10! * (20-10)! ) 184756 > > if (FALSE) + {" + There are 184756 possible configurations placing 10 of the 20 observations into Sample 1 and the other + 10 observations going into Sample 2. Each configuration is equally likely to occur under the null hypothesis. + The logic of the permutation test is quite similar to the logic of the t-test, but no assumption of normality + is needed. If the sampled case results in a metric (t statistic) that is in the most extreme 5% of possible + results, then we reject the null hypothesis of no difference between the means. + The advantage of the permutation test is it does not make any assumption about normality, and in fact, doesn't + make any assumption at all about parent distributions. The disadvantage of a permutation test is the + number of permutations that must be generated. One possible solution is if the number of permutations is + large then take a random sample (without replacement) of possible permuations + "} > > set.seed(13254) #for reproducibility of results > > R scores t.values > for (i in 1:R) { > index #we took a random sample of the vector 1 to 20, without replacement, programming is required to assure no > # and stored that in an object called index. These index values were then used > # to extract the first resampled group of scores. duplicates of the possible configurations. > #The second group of scores was extracted using a trick, which in words goes, > #for group two take all the scores that are not indexed by the values in index. > # In other words, put the remaining scores in group two. Then we did the t-test and stored > # the test statistic (t-value). It now remains to compare those simulated t-values to the one > # we obtained above when we did the actual t-test. > # > t.values = abs(t.values) # for a two-tailed test > #compute the proportion of abs(t-values) >= 0.048008 > #note: 0.048008 was the computed t-statistic from the two-sample test on our data > mean(t.values >=0.048008) #a mean of 0s and 1s is equal to the number of 1s 0.9420317 > # > # compare to the p-value from the two-sample t-test that was equal to 0.9622 > #The resampling scheme tells us that a little more than 0.9420 of the simulations gave us a abs(t.value) > #as extreme or more extreme than the abs(t* observed in the experiement). > #This is the p-value resulting from the permutation test. > #We can see that it is pretty close to the p-value obtained in the actual t-test above (0.9622). > ##-------------------------------------------------------------## > dev.off() null device 1 >