###################################################################### # Computer notes Statistics for Laboratory Scientists (2) # Lecture 3 Johns Hopkins University ###################################################################### # These notes provide R code related to the course lectures, to # further illustrate the topics in the lectures and to assist the # students to learn R. # # Lines beginning with the symbol '#' are comments in R. All other # lines contain code. # # In R for Windows, you may wish to open this file from the menu bar # (File:Display file); you can then easily copy commands into the # command window. (Use the mouse to highlight one or more lines; then # right-click and select "Paste to console".) ###################################################################### ###################################################################### # example 1 ###################################################################### x <- rbind(c(18, 2), c(11, 9)) # marginal totals rs <- apply(x, 1, sum) # sum of each row cs <- apply(x, 2, sum) # sum of each column n <- sum(x) # expected counts expected <- outer(rs,cs,"*")/n # chi-square statistic chi <- sum((x-expected)^2/expected) # result = 6.14 # P-value 1-pchisq(chi,1) # value=0.013 # lrt statistic lrt <- 2*sum(x * log(x/expected)) # result = 6.52 # P-value 1-pchisq(lrt,1)# value=0.011 # Fisher's exact test fisher.test(x) # p-value = 3.1% ###################################################################### # example 2 ###################################################################### x <- rbind(c(9,9), c(20, 62)) # marginal totals rs <- apply(x, 1, sum) # sum of each row cs <- apply(x, 2, sum) # sum of each column n <- sum(x) # expected counts expected <- outer(rs,cs,"*")/n # chi-square statistic chi <- sum((x-expected)^2/expected) # result = 4.70 # P-value 1-pchisq(chi,1) # value=0.030 # lrt statistic lrt <- 2*sum(x * log(x/expected)) # result = 4.37 # P-value 1-pchisq(lrt,1)# value=0.037 # Fisher's exact test fisher.test(x) # p-value = 4.4% ###################################################################### # Hypergeometric distribution ###################################################################### # Consider an urn with 30 balls, of which 12 are white and 18 are black # Suppose we draw 10 balls *without* replacement and count the number # of white balls obtained. ###################################################################### # Simulate this rhyper(1, 12, 18, 10) # Simulate this 100 times rhyper(100, 12, 18, 10) # Probability of getting no white balls in the sample dhyper(0, 12, 18, 10) # Probability of getting exactly 2 white balls in the sample dhyper(2, 12, 18, 10) # Probability of getting 2 or fewer white balls in the sample phyper(2, 12, 18, 10) # 2.5th and 97.5th %iles of number of white balls drawn qhyper(c(0.025, 0.975), 12, 18, 10) ###################################################################### # Paired data example ###################################################################### # The data z <- rbind(c(9,9),c(20,62)) # McNemar's test chi <- (z[1,2]-z[2,1])^2 / (z[1,2]+z[2,1]) # value = 4.17 1 - pchisq(chi, 1) # P = 4.1% # An exact test binom.test(z[1,2], z[1,2]+z[2,1]) # P = 6.1% ################## # End of comp03.R ##################