Oconnor_Barbra
Instructions
Submit both the .Rmd and .html files for grading.
Please delete the Instructions shown above prior to submitting your .Rmd and .html files.
Test Items (50 points total)
(1) R has probability functions available for use (Kabacoff, Section 5.2.3). Using one distribution to approximate another is not uncommon.
(1)(a) (4 points) The Poisson distribution may be used to approximate the binomial distribution if n > 20 and np < 7. Estimate the following binomial probabilities using dpois() and ppois() with probability p = 0.05, and n = 100. Then, estimate the same probabilities using dbinom() and pbinom(). Show the numerical results of your calculations.
- The probability of exactly 0 successes.
dpois(0,lambda = 5)## [1] 0.006737947ppois(0,lambda = 5)## [1] 0.006737947dbinom(0,100,0.05)## [1] 0.005920529pbinom(0,100,0.05)## [1] 0.005920529- The probability of fewer than 6 successes.
sum(dpois(0:5,lambda = 5))## [1] 0.6159607ppois(5,lambda = 5)## [1] 0.6159607sum(dbinom(0:5,100,0.05))## [1] 0.6159991pbinom(5,100,0.05)## [1] 0.6159991(1)(b) (2 points) Generate side-by-side barplots using par(mfrow = c(1,2)) or grid.arrange(). The left barplot will show Poisson probabilties for outcomes ranging from 0 to 10. The right barplot will show binomial probabilities for outcomes ranging from 0 to 10. Use p = 0.05 and n = 100. Title each plot, present in color and assign names to the bar; i.e. x-axis value labels.
dp<-(dpois(0:10,lambda = 5))
dp## [1] 0.006737947 0.033689735 0.084224337 0.140373896 0.175467370
## [6] 0.175467370 0.146222808 0.104444863 0.065278039 0.036265577
## [11] 0.018132789db<-(dbinom(0:10,100,0.05))
db## [1] 0.005920529 0.031160680 0.081181772 0.139575678 0.178142642
## [6] 0.180017827 0.150014856 0.106025537 0.064870888 0.034901296
## [11] 0.016715884par(mfrow=c(1,2))
barplot(dp,main ="Poisson Probability",col="aquamarine3",xlab = "Outcome",ylab ="Probability")
barplot(db,main = "Binomial Probability",col = "coral",xlab = "Outcome",ylab = "Probability")
(1)(c) For this problem refer to Sections 5.2 of Business Statistics. A discrete random variable has outcomes: 0, 1, 2, 3, 4, 5, 6. The corresponding probabilities in sequence with the outcomes are: 0.215, 0.230, 0.240, 0.182, 0.130, 0.003, 0.001. In other words, the probabilty of obtaining “0” is 0.215.
- (2 points) Calculate the expected value and variance for this distribution using the general formula for mean and variance of a discrete distribution. To do this, you will need to use integer values from 0 to 6 as outcomes along with the corresponding probabilities. Round your answer to 2 decimal places.
x<-0:6
prob<-c(0.215,0.230,0.240,0.182,0.130,0.003,0.001)
meanx<-round(sum(x*prob),digits=2)
meanx## [1] 1.8varx<-round(((sum((x^2)*prob))-(meanx^2)),digits=2)
varx## [1] 1.78- (2 points) Use the cumsum() function and plot the cumulative probabilties versus the corresponding outcomes. Detemine the value of the median for this distribution and show on this plot.
cumprob<-cumsum(prob)
cumprob## [1] 0.215 0.445 0.685 0.867 0.997 1.000 1.001plot(x,cumprob,pch=22,main = "Cum. Prob. vs Corr. Outcomes",xlab="outcomes",ylab="Cum prob")
text(1.8,0.56,"The Median is assumed to be 1.8")
#mean assumed to be (1.8,0) based on the answer to Prob. #1.c(i) above(2) Conditional probabilities appear in many contexts and in particular are used by Bayes’ Theorem. Correlations are another means for evaluating dependency between variables. The dataset “faithful”" is part of the “datasets” package and may be loaded with the statement data(faithful). It contains 272 observations of 2 variables; waiting time between eruptions (in minutes) and the duration of the eruption (in minutes) for the Old Faithful geyser in Yellowstone National Park.
(2)(a) (2 points) Load the “faithful” and present summary statistics and a histogram of waiting times. Additionally, compute the empirical conditional probability of an eruption less than 3.0 minutes, if the waiting time exceeds 70 minutes.
data("faithful")
summary(faithful)## eruptions waiting
## Min. :1.600 Min. :43.0
## 1st Qu.:2.163 1st Qu.:58.0
## Median :4.000 Median :76.0
## Mean :3.488 Mean :70.9
## 3rd Qu.:4.454 3rd Qu.:82.0
## Max. :5.100 Max. :96.0str(faithful)## 'data.frame': 272 obs. of 2 variables:
## $ eruptions: num 3.6 1.8 3.33 2.28 4.53 ...
## $ waiting : num 79 54 74 62 85 55 88 85 51 85 ...hist(faithful$waiting)
faithful## eruptions waiting
## 1 3.600 79
## 2 1.800 54
## 3 3.333 74
## 4 2.283 62
## 5 4.533 85
## 6 2.883 55
## 7 4.700 88
## 8 3.600 85
## 9 1.950 51
## 10 4.350 85
## 11 1.833 54
## 12 3.917 84
## 13 4.200 78
## 14 1.750 47
## 15 4.700 83
## 16 2.167 52
## 17 1.750 62
## 18 4.800 84
## 19 1.600 52
## 20 4.250 79
## 21 1.800 51
## 22 1.750 47
## 23 3.450 78
## 24 3.067 69
## 25 4.533 74
## 26 3.600 83
## 27 1.967 55
## 28 4.083 76
## 29 3.850 78
## 30 4.433 79
## 31 4.300 73
## 32 4.467 77
## 33 3.367 66
## 34 4.033 80
## 35 3.833 74
## 36 2.017 52
## 37 1.867 48
## 38 4.833 80
## 39 1.833 59
## 40 4.783 90
## 41 4.350 80
## 42 1.883 58
## 43 4.567 84
## 44 1.750 58
## 45 4.533 73
## 46 3.317 83
## 47 3.833 64
## 48 2.100 53
## 49 4.633 82
## 50 2.000 59
## 51 4.800 75
## 52 4.716 90
## 53 1.833 54
## 54 4.833 80
## 55 1.733 54
## 56 4.883 83
## 57 3.717 71
## 58 1.667 64
## 59 4.567 77
## 60 4.317 81
## 61 2.233 59
## 62 4.500 84
## 63 1.750 48
## 64 4.800 82
## 65 1.817 60
## 66 4.400 92
## 67 4.167 78
## 68 4.700 78
## 69 2.067 65
## 70 4.700 73
## 71 4.033 82
## 72 1.967 56
## 73 4.500 79
## 74 4.000 71
## 75 1.983 62
## 76 5.067 76
## 77 2.017 60
## 78 4.567 78
## 79 3.883 76
## 80 3.600 83
## 81 4.133 75
## 82 4.333 82
## 83 4.100 70
## 84 2.633 65
## 85 4.067 73
## 86 4.933 88
## 87 3.950 76
## 88 4.517 80
## 89 2.167 48
## 90 4.000 86
## 91 2.200 60
## 92 4.333 90
## 93 1.867 50
## 94 4.817 78
## 95 1.833 63
## 96 4.300 72
## 97 4.667 84
## 98 3.750 75
## 99 1.867 51
## 100 4.900 82
## 101 2.483 62
## 102 4.367 88
## 103 2.100 49
## 104 4.500 83
## 105 4.050 81
## 106 1.867 47
## 107 4.700 84
## 108 1.783 52
## 109 4.850 86
## 110 3.683 81
## 111 4.733 75
## 112 2.300 59
## 113 4.900 89
## 114 4.417 79
## 115 1.700 59
## 116 4.633 81
## 117 2.317 50
## 118 4.600 85
## 119 1.817 59
## 120 4.417 87
## 121 2.617 53
## 122 4.067 69
## 123 4.250 77
## 124 1.967 56
## 125 4.600 88
## 126 3.767 81
## 127 1.917 45
## 128 4.500 82
## 129 2.267 55
## 130 4.650 90
## 131 1.867 45
## 132 4.167 83
## 133 2.800 56
## 134 4.333 89
## 135 1.833 46
## 136 4.383 82
## 137 1.883 51
## 138 4.933 86
## 139 2.033 53
## 140 3.733 79
## 141 4.233 81
## 142 2.233 60
## 143 4.533 82
## 144 4.817 77
## 145 4.333 76
## 146 1.983 59
## 147 4.633 80
## 148 2.017 49
## 149 5.100 96
## 150 1.800 53
## 151 5.033 77
## 152 4.000 77
## 153 2.400 65
## 154 4.600 81
## 155 3.567 71
## 156 4.000 70
## 157 4.500 81
## 158 4.083 93
## 159 1.800 53
## 160 3.967 89
## 161 2.200 45
## 162 4.150 86
## 163 2.000 58
## 164 3.833 78
## 165 3.500 66
## 166 4.583 76
## 167 2.367 63
## 168 5.000 88
## 169 1.933 52
## 170 4.617 93
## 171 1.917 49
## 172 2.083 57
## 173 4.583 77
## 174 3.333 68
## 175 4.167 81
## 176 4.333 81
## 177 4.500 73
## 178 2.417 50
## 179 4.000 85
## 180 4.167 74
## 181 1.883 55
## 182 4.583 77
## 183 4.250 83
## 184 3.767 83
## 185 2.033 51
## 186 4.433 78
## 187 4.083 84
## 188 1.833 46
## 189 4.417 83
## 190 2.183 55
## 191 4.800 81
## 192 1.833 57
## 193 4.800 76
## 194 4.100 84
## 195 3.966 77
## 196 4.233 81
## 197 3.500 87
## 198 4.366 77
## 199 2.250 51
## 200 4.667 78
## 201 2.100 60
## 202 4.350 82
## 203 4.133 91
## 204 1.867 53
## 205 4.600 78
## 206 1.783 46
## 207 4.367 77
## 208 3.850 84
## 209 1.933 49
## 210 4.500 83
## 211 2.383 71
## 212 4.700 80
## 213 1.867 49
## 214 3.833 75
## 215 3.417 64
## 216 4.233 76
## 217 2.400 53
## 218 4.800 94
## 219 2.000 55
## 220 4.150 76
## 221 1.867 50
## 222 4.267 82
## 223 1.750 54
## 224 4.483 75
## 225 4.000 78
## 226 4.117 79
## 227 4.083 78
## 228 4.267 78
## 229 3.917 70
## 230 4.550 79
## 231 4.083 70
## 232 2.417 54
## 233 4.183 86
## 234 2.217 50
## 235 4.450 90
## 236 1.883 54
## 237 1.850 54
## 238 4.283 77
## 239 3.950 79
## 240 2.333 64
## 241 4.150 75
## 242 2.350 47
## 243 4.933 86
## 244 2.900 63
## 245 4.583 85
## 246 3.833 82
## 247 2.083 57
## 248 4.367 82
## 249 2.133 67
## 250 4.350 74
## 251 2.200 54
## 252 4.450 83
## 253 3.567 73
## 254 4.500 73
## 255 4.150 88
## 256 3.817 80
## 257 3.917 71
## 258 4.450 83
## 259 2.000 56
## 260 4.283 79
## 261 4.767 78
## 262 4.533 84
## 263 1.850 58
## 264 4.250 83
## 265 1.983 43
## 266 2.250 60
## 267 4.750 75
## 268 4.117 81
## 269 2.150 46
## 270 4.417 90
## 271 1.817 46
## 272 4.467 74Wait70<-subset(faithful,waiting>70)
Wait70## eruptions waiting
## 1 3.600 79
## 3 3.333 74
## 5 4.533 85
## 7 4.700 88
## 8 3.600 85
## 10 4.350 85
## 12 3.917 84
## 13 4.200 78
## 15 4.700 83
## 18 4.800 84
## 20 4.250 79
## 23 3.450 78
## 25 4.533 74
## 26 3.600 83
## 28 4.083 76
## 29 3.850 78
## 30 4.433 79
## 31 4.300 73
## 32 4.467 77
## 34 4.033 80
## 35 3.833 74
## 38 4.833 80
## 40 4.783 90
## 41 4.350 80
## 43 4.567 84
## 45 4.533 73
## 46 3.317 83
## 49 4.633 82
## 51 4.800 75
## 52 4.716 90
## 54 4.833 80
## 56 4.883 83
## 57 3.717 71
## 59 4.567 77
## 60 4.317 81
## 62 4.500 84
## 64 4.800 82
## 66 4.400 92
## 67 4.167 78
## 68 4.700 78
## 70 4.700 73
## 71 4.033 82
## 73 4.500 79
## 74 4.000 71
## 76 5.067 76
## 78 4.567 78
## 79 3.883 76
## 80 3.600 83
## 81 4.133 75
## 82 4.333 82
## 85 4.067 73
## 86 4.933 88
## 87 3.950 76
## 88 4.517 80
## 90 4.000 86
## 92 4.333 90
## 94 4.817 78
## 96 4.300 72
## 97 4.667 84
## 98 3.750 75
## 100 4.900 82
## 102 4.367 88
## 104 4.500 83
## 105 4.050 81
## 107 4.700 84
## 109 4.850 86
## 110 3.683 81
## 111 4.733 75
## 113 4.900 89
## 114 4.417 79
## 116 4.633 81
## 118 4.600 85
## 120 4.417 87
## 123 4.250 77
## 125 4.600 88
## 126 3.767 81
## 128 4.500 82
## 130 4.650 90
## 132 4.167 83
## 134 4.333 89
## 136 4.383 82
## 138 4.933 86
## 140 3.733 79
## 141 4.233 81
## 143 4.533 82
## 144 4.817 77
## 145 4.333 76
## 147 4.633 80
## 149 5.100 96
## 151 5.033 77
## 152 4.000 77
## 154 4.600 81
## 155 3.567 71
## 157 4.500 81
## 158 4.083 93
## 160 3.967 89
## 162 4.150 86
## 164 3.833 78
## 166 4.583 76
## 168 5.000 88
## 170 4.617 93
## 173 4.583 77
## 175 4.167 81
## 176 4.333 81
## 177 4.500 73
## 179 4.000 85
## 180 4.167 74
## 182 4.583 77
## 183 4.250 83
## 184 3.767 83
## 186 4.433 78
## 187 4.083 84
## 189 4.417 83
## 191 4.800 81
## 193 4.800 76
## 194 4.100 84
## 195 3.966 77
## 196 4.233 81
## 197 3.500 87
## 198 4.366 77
## 200 4.667 78
## 202 4.350 82
## 203 4.133 91
## 205 4.600 78
## 207 4.367 77
## 208 3.850 84
## 210 4.500 83
## 211 2.383 71
## 212 4.700 80
## 214 3.833 75
## 216 4.233 76
## 218 4.800 94
## 220 4.150 76
## 222 4.267 82
## 224 4.483 75
## 225 4.000 78
## 226 4.117 79
## 227 4.083 78
## 228 4.267 78
## 230 4.550 79
## 233 4.183 86
## 235 4.450 90
## 238 4.283 77
## 239 3.950 79
## 241 4.150 75
## 243 4.933 86
## 245 4.583 85
## 246 3.833 82
## 248 4.367 82
## 250 4.350 74
## 252 4.450 83
## 253 3.567 73
## 254 4.500 73
## 255 4.150 88
## 256 3.817 80
## 257 3.917 71
## 258 4.450 83
## 260 4.283 79
## 261 4.767 78
## 262 4.533 84
## 264 4.250 83
## 267 4.750 75
## 268 4.117 81
## 270 4.417 90
## 272 4.467 74erupt3<-subset(Wait70,eruptions<3)
erupt3## eruptions waiting
## 211 2.383 71condprob<-nrow(erupt3)/nrow(Wait70)
condprob## [1] 0.006060606- (2 points) Identify any observations in “faithful” for which the waiting time exceeds 70 minutes and the eruptions are less than 3.0 minutes. List and show any such observations in a distinct color on a scatterplot of all eruption (vertical axis) and waiting times (horizontal axis). Include a horizontal line at eruption = 3.0, and a vertical line at waiting time = 70. Add a title and appropriate text.
waiting = faithful$waiting
duration = faithful$eruptions
plot(waiting, duration,main="Faithful Observations", xlab = "Time Waited",ylab = "Eruption Duration")
abline(h=3,v=70,col="blue")
- (1 point) What does the plot suggest about the relationship between eruption time and waiting time?
Answer: (There is generally a positive, linear relationship between eruption time and waiting time with Wait Times Increasing when the Duration of the Eruption also increases )
(2)(b) (3 points) Past research indicates that the waiting times between consecutive eruptions are not independent. This problem will check to see if there is evidence of this. Form consecutive pairs of waiting times. In other words, pair the first and second waiting times, pair the third and fourth waiting times, and so forth. There are 136 resulting consecutive pairs of waiting times. Form a data frame with the first column containing the first waiting time in a pair and the second column with the second waiting time in a pair. Plot the pairs with the second member of a pair on the vertical axis and the first member on the horizontal axis.
One way to do this is to pass the vector of waiting times - faithful$waiting - to matrix(), specifying 2 columns for our matrix, with values organized by row; i.e. byrow = TRUE.
dim(faithful)## [1] 272 2colnames(faithful)## [1] "eruptions" "waiting"mwait<-matrix(waiting,ncol=2,byrow=TRUE)
mwait## [,1] [,2]
## [1,] 79 54
## [2,] 74 62
## [3,] 85 55
## [4,] 88 85
## [5,] 51 85
## [6,] 54 84
## [7,] 78 47
## [8,] 83 52
## [9,] 62 84
## [10,] 52 79
## [11,] 51 47
## [12,] 78 69
## [13,] 74 83
## [14,] 55 76
## [15,] 78 79
## [16,] 73 77
## [17,] 66 80
## [18,] 74 52
## [19,] 48 80
## [20,] 59 90
## [21,] 80 58
## [22,] 84 58
## [23,] 73 83
## [24,] 64 53
## [25,] 82 59
## [26,] 75 90
## [27,] 54 80
## [28,] 54 83
## [29,] 71 64
## [30,] 77 81
## [31,] 59 84
## [32,] 48 82
## [33,] 60 92
## [34,] 78 78
## [35,] 65 73
## [36,] 82 56
## [37,] 79 71
## [38,] 62 76
## [39,] 60 78
## [40,] 76 83
## [41,] 75 82
## [42,] 70 65
## [43,] 73 88
## [44,] 76 80
## [45,] 48 86
## [46,] 60 90
## [47,] 50 78
## [48,] 63 72
## [49,] 84 75
## [50,] 51 82
## [51,] 62 88
## [52,] 49 83
## [53,] 81 47
## [54,] 84 52
## [55,] 86 81
## [56,] 75 59
## [57,] 89 79
## [58,] 59 81
## [59,] 50 85
## [60,] 59 87
## [61,] 53 69
## [62,] 77 56
## [63,] 88 81
## [64,] 45 82
## [65,] 55 90
## [66,] 45 83
## [67,] 56 89
## [68,] 46 82
## [69,] 51 86
## [70,] 53 79
## [71,] 81 60
## [72,] 82 77
## [73,] 76 59
## [74,] 80 49
## [75,] 96 53
## [76,] 77 77
## [77,] 65 81
## [78,] 71 70
## [79,] 81 93
## [80,] 53 89
## [81,] 45 86
## [82,] 58 78
## [83,] 66 76
## [84,] 63 88
## [85,] 52 93
## [86,] 49 57
## [87,] 77 68
## [88,] 81 81
## [89,] 73 50
## [90,] 85 74
## [91,] 55 77
## [92,] 83 83
## [93,] 51 78
## [94,] 84 46
## [95,] 83 55
## [96,] 81 57
## [97,] 76 84
## [98,] 77 81
## [99,] 87 77
## [100,] 51 78
## [101,] 60 82
## [102,] 91 53
## [103,] 78 46
## [104,] 77 84
## [105,] 49 83
## [106,] 71 80
## [107,] 49 75
## [108,] 64 76
## [109,] 53 94
## [110,] 55 76
## [111,] 50 82
## [112,] 54 75
## [113,] 78 79
## [114,] 78 78
## [115,] 70 79
## [116,] 70 54
## [117,] 86 50
## [118,] 90 54
## [119,] 54 77
## [120,] 79 64
## [121,] 75 47
## [122,] 86 63
## [123,] 85 82
## [124,] 57 82
## [125,] 67 74
## [126,] 54 83
## [127,] 73 73
## [128,] 88 80
## [129,] 71 83
## [130,] 56 79
## [131,] 78 84
## [132,] 58 83
## [133,] 43 60
## [134,] 75 81
## [135,] 46 90
## [136,] 46 74plot(y=mwait[,2],x=mwait[,1],xlab="First Column",ylab = "Second Column")
(2)(c) (2) Test the hypothesis of independence with a two-sided test at the 5% level using the Kendall correlation coefficient.
library(Kendall)
cor.test(mwait[,1],mwait[,2],alternative = "two.sided",method = "kendall", conf.level = 0.95)##
## Kendall's rank correlation tau
##
## data: mwait[, 1] and mwait[, 2]
## z = -4.9482, p-value = 7.489e-07
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## -0.2935579#Null hypothesis = data are not independent; we reject the Null hypothesis(3) Performing hypothesis tests using random samples is fundamental to statistical inference. The first part of this problem involves comparing two different diets. Using “ChickWeight” data available in the base R, “datasets” package, execute the following code to prepare a data frame for analysis.
# load "ChickWeight" dataset
data(ChickWeight)
# Create T | F vector indicating observations with Time == 21 and Diet == "1" OR "3"
index <- ChickWeight$Time == 21 & (ChickWeight$Diet == "1" | ChickWeight$Diet == "3")
# Create data frame, "result," with the weight and Diet of those observations with "TRUE" "index"" values
result <- subset(ChickWeight[index, ], select = c(weight, Diet))
# Encode "Diet" as a factor
result$Diet <- factor(result$Diet)
str(result) ## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 26 obs. of 2 variables:
## $ weight: num 205 215 202 157 223 157 305 98 124 175 ...
## $ Diet : Factor w/ 2 levels "1","3": 1 1 1 1 1 1 1 1 1 1 ...The data frame, “result”, has chick weights for two diets, identified as diet “1” and “3”. Use the data frame, “result,” to complete the following item.
(3)(a) (2 points) Display two side-by-side vertical boxplots using par(mfrow = c(1,2)). One boxplot would display diet “1” and the other diet “3”.
par(mfrow=c(1,2))
result## weight Diet
## 12 205 1
## 24 215 1
## 36 202 1
## 48 157 1
## 60 223 1
## 72 157 1
## 84 305 1
## 107 98 1
## 119 124 1
## 131 175 1
## 143 205 1
## 155 96 1
## 167 266 1
## 194 142 1
## 208 157 1
## 220 117 1
## 352 256 3
## 364 305 3
## 376 147 3
## 388 341 3
## 400 373 3
## 412 220 3
## 424 178 3
## 436 290 3
## 448 272 3
## 460 321 3result1<-subset(result,Diet == "1")
result1## weight Diet
## 12 205 1
## 24 215 1
## 36 202 1
## 48 157 1
## 60 223 1
## 72 157 1
## 84 305 1
## 107 98 1
## 119 124 1
## 131 175 1
## 143 205 1
## 155 96 1
## 167 266 1
## 194 142 1
## 208 157 1
## 220 117 1boxplot(result1$weight)
result3<-subset(result,Diet == "3")
boxplot(result3$weight)
(3)(b) (2 points) Use the “weight” data for the two diets to test the null hypothesis of equal population mean weights for the two diets. Test at the 95% confidence level with a two-sided t-test. This can be done using t.test() in R. Assume equal variances. Display the results of t.test().
t.test(result1$weight,result3$weight,alternative = "two.sided",conf.level = 0.95)##
## Welch Two Sample t-test
##
## data: result1$weight and result3$weight
## t = -3.4293, df = 16.408, p-value = 0.003337
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -149.64644 -35.45356
## sample estimates:
## mean of x mean of y
## 177.75 270.30Working with paired data is another common statistical activity. The “ChickWeight” data will be used to illustrate how the weight gain from day 20 to 21 may be analyzed. Use the following code to prepare pre- and post-data from Diet == “3” for analysis.
# load "ChickWeight" dataset
data(ChickWeight)
# Create T | F vector indicating observations with Diet == "3"
index <- ChickWeight$Diet == "3"
# Create vector of "weight" for observations where Diet == "3" and Time == 20
pre <- subset(ChickWeight[index, ], Time == 20, select = weight)$weight
# Create vector of "weight" for observations where Diet == "3" and Time == 21
post <- subset(ChickWeight[index, ], Time == 21, select = weight)$weight
# The pre and post values are paired, each pair corresponding to an individual chick.
cbind(pre, post)## pre post
## [1,] 235 256
## [2,] 291 305
## [3,] 156 147
## [4,] 327 341
## [5,] 361 373
## [6,] 225 220
## [7,] 169 178
## [8,] 280 290
## [9,] 250 272
## [10,] 295 321(3)(c) (2 points) Present a scatterplot of the variable “post” as a function of the variable “pre”. Include a diagonal line with zero intercept and slope equal to one. Title and label the variables in this scatterplot.
plot(pre,post,col="blue",xlab="Variable Pre",ylab="Variable Post",
main="Scatterplot Pre~Post",pch=15)
abline(0,1,col="red")
(3)(d) (4 points) Calculate and present a one-sided, 95% confidence interval for the average weight gain from day 20 to day 21. Write the code for the paired t-test and for determination of the confidence interval endpoints. **Do not use *t.test()**, although you may check your answers using this function. Present the resulting test statistic value, critical value, p-value and confidence interval.
mean_pre<-mean(pre)
mean_pre## [1] 258.9mean_post<-mean(post)
mean_post## [1] 270.3mean_weight_gain<-mean_post-mean_pre
mean_weight_gain## [1] 11.4t.test(pre,post,alternate="two.sided",conf.level = 0.95)##
## Welch Two Sample t-test
##
## data: pre and post
## t = -0.37209, df = 17.846, p-value = 0.7142
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -75.8069 53.0069
## sample estimates:
## mean of x mean of y
## 258.9 270.3n<-length(pre)
n## [1] 10df<-n-1
df## [1] 9t_lower<-qt(0.025,df)
t_lower## [1] -2.262157t_upper<-qt(0.975,df)
t_upper## [1] 2.262157std.dev<-sqrt(sum(mean_weight_gain^2)/(n-1))
std.dev## [1] 3.8std_error<-std.dev/sqrt(n)
std_error## [1] 1.201666lowerci<-mean_weight_gain+(std_error*t_lower)
lowerci## [1] 8.681644upperci<-mean_weight_gain+(std_error*t_upper)
upperci## [1] 14.11836(4) Statistical inference depends on using a sampling distribution for a statistic in order to make confidence statements about unknown population parameters. The Central Limit Theorem is used to justify use of the normal distribution as a sampling distribution for statistical inference. Using Nile River flow data from 1871 to 1970, this problem demonstrates sampling distribution convergence to normality. Use the code below to prepare the data. Refer to this example when completing (4)(c) below.
data(Nile)
m <- mean(Nile)
std <- sd(Nile)
x <- seq(from = 400, to = 1400, by = 1)
hist(Nile, freq = FALSE, col = "darkblue", xlab = "Flow",
main = "Histogram of Nile River Flows, 1871 to 1970")
curve(dnorm(x, mean = m, sd = std), col = "orange", lwd = 2, add = TRUE)
(4)(a) (2 points) Using Nile River flow data and the “moments” package, calculate skewness and kurtosis. Present a QQ plot and boxplot of the flow data side-by-side using qqnorm(), qqline() and boxplot(); par(mfrow = c(1, 2)) may be used to locate the plots side-by-side. Add features to these displays as you choose.
library(moments)
skewness1<-skewness(Nile)
skewness1## [1] 0.3223697kurtosis1<-kurtosis(Nile)
kurtosis1## [1] 2.695093par(mfrow=c(1,2))
qqnorm(Nile,main = "QQ PLOT")
qqline(Nile)
boxplot(Nile,main="Boxplot")
(4)(b) (4 points) Using set.seed(124) and the Nile data, generate 1000 random samples of size n = 16, with replacement. For each sample drawn, calculate and store the sample mean. This can be done with a for-loop and use of the sample() function. Label the resulting 1000 mean values as “sample1”. Repeat these steps using set.seed(127) - a different “seed” - and samples of size n = 64. Label these 1000 mean values as “sample2”. Compute and present the means, sample standard deviations and sample variances for “sample1” and “sample2” in a table with the first row for “sample1”, the second row for “sample2” and the columns labled for each statistic.
set.seed(124)
sample1<-rep(1:1000,0)
for(i in 1:1000) {
sample1[i]<-mean(sample(Nile,16,replace=TRUE))
}
set.seed(127)
sample2<-rep(1:1000,0)
for(i in 1:1000) {
sample2[i]<-mean(sample(Nile,64,replace = TRUE))
}
row_names<-c("sample1","sample2")
col_names<-c("mean","sample std dev","sample variance")
matrix(c(mean(sample1),mean(sample2),sd(sample1),sd(sample2),var(sample1),var(sample2)),nrow = 2,ncol = 3,dimnames = list(row_names,col_names))## mean sample std dev sample variance
## sample1 918.7349 41.57159 1728.197
## sample2 919.4785 20.07329 402.937(4)(c) (4 points) Present side-by-side histograms of “sample1” and “sample2” with the normal density curve superimposed. To prepare comparable histograms, it will be necessary to use “freq = FALSE” and to maintain the same x-axis with “xlim = c(750, 1050)”, and the same y-axis with “ylim = c(0, 0.025).” To superimpose separate density functions, you will need to use the mean and standard deviation for each “sample” - each histogram - separately.
par(mfrow=c(1,2))
hist(sample1,freq = FALSE,xlim = c(750,1040),ylim = c(0,0.025))
curve(dnorm(x,mean=mean(sample1),sd=sd(sample1)),add=TRUE,col="red",lty=2,lwd=2)
hist(sample2,freq = FALSE,xlim = c(750,1040),ylim = c(0,0.025))
curve(dnorm(x,mean=mean(sample2),sd=sd(sample2)),add=TRUE,col="red",lty=2,lwd=2)
(5) This problem deals with contingency table analysis. This is an example of categorical data analysis (see Kabacoff, pp. 145-151). The “warpbreaks” dataset gives the number of warp breaks per loom, where a loom corresponds to a fixed length of yarn. There are 54 observations on 3 variables: breaks (numeric, the number of breaks), wool (factor, type of wool: A or B), and tension (factor, low L, medium M and high H). These data have been studied and used for example elsewhere. For the purposes of this problem, we will focus on the relationship between breaks and tension using contingency table analysis.
(5)(a)(3 points) warpbreaks is part of the “datasets” package and may be loaded via data(warpbreaks). Load “warpbreaks” and present the structure using str(). Calculate the median number of breaks for the entire dataset, disregarding “tension” and “wool”. Define this median value as “median_breaks”. Present a histogram of the number of breaks with the location of the median indicated.
Create a new variable “number” as follows: for each value of “breaks”, classify the number of breaks as either strictly below “median_breaks”, or the alternative. Convert the “above”|“below” classifications to a factor, and combine with the dataset warpbreaks. Present a summary of the augmented dataset using summary(). Present a contingency table of the frequency of breaks using the two variables “tension” and “number”. There should be six cells in this table.
data("warpbreaks")
str(warpbreaks)## 'data.frame': 54 obs. of 3 variables:
## $ breaks : num 26 30 54 25 70 52 51 26 67 18 ...
## $ wool : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
## $ tension: Factor w/ 3 levels "L","M","H": 1 1 1 1 1 1 1 1 1 2 ...median_breaks<-median(warpbreaks$breaks)
median_breaks## [1] 26hist(warpbreaks$breaks)
abline(v = median(warpbreaks$breaks),
col = "royalblue",
lwd = 2)
number<-ifelse(warpbreaks$breaks<median_breaks,"below","above")
number## [1] "above" "above" "above" "below" "above" "above" "above" "above"
## [9] "above" "below" "below" "above" "below" "below" "below" "above"
## [17] "above" "above" "above" "below" "below" "below" "below" "above"
## [25] "above" "below" "above" "above" "below" "above" "below" "above"
## [33] "above" "above" "below" "above" "above" "above" "below" "below"
## [41] "above" "above" "below" "above" "above" "below" "below" "below"
## [49] "below" "below" "below" "below" "below" "above"warpbreaks2<-cbind(warpbreaks,number)
warpbreaks2## breaks wool tension number
## 1 26 A L above
## 2 30 A L above
## 3 54 A L above
## 4 25 A L below
## 5 70 A L above
## 6 52 A L above
## 7 51 A L above
## 8 26 A L above
## 9 67 A L above
## 10 18 A M below
## 11 21 A M below
## 12 29 A M above
## 13 17 A M below
## 14 12 A M below
## 15 18 A M below
## 16 35 A M above
## 17 30 A M above
## 18 36 A M above
## 19 36 A H above
## 20 21 A H below
## 21 24 A H below
## 22 18 A H below
## 23 10 A H below
## 24 43 A H above
## 25 28 A H above
## 26 15 A H below
## 27 26 A H above
## 28 27 B L above
## 29 14 B L below
## 30 29 B L above
## 31 19 B L below
## 32 29 B L above
## 33 31 B L above
## 34 41 B L above
## 35 20 B L below
## 36 44 B L above
## 37 42 B M above
## 38 26 B M above
## 39 19 B M below
## 40 16 B M below
## 41 39 B M above
## 42 28 B M above
## 43 21 B M below
## 44 39 B M above
## 45 29 B M above
## 46 20 B H below
## 47 21 B H below
## 48 24 B H below
## 49 17 B H below
## 50 13 B H below
## 51 15 B H below
## 52 15 B H below
## 53 16 B H below
## 54 28 B H abovesummary(warpbreaks2)## breaks wool tension number
## Min. :10.00 A:27 L:18 above:29
## 1st Qu.:18.25 B:27 M:18 below:25
## Median :26.00 H:18
## Mean :28.15
## 3rd Qu.:34.00
## Max. :70.00contingency<-table(warpbreaks2$tension,warpbreaks2$number)
contingency##
## above below
## L 14 4
## M 10 8
## H 5 13(5)(b)(2 points) Using the table constructed in (5)(a), test at the 5% level the null hypothesis of independence using the uncorrected chisq.test() (Black, Business Statistics, Section 16.2). Show the results of this test and state your conclusions.
chisq.test(contingency)##
## Pearson's Chi-squared test
##
## data: contingency
## X-squared = 9.0869, df = 2, p-value = 0.01064#As the resulting p-value is small, reject the Null hypothesis and the variables are not independent(5)(c) (5 points) Write a function that computes the uncorrected Pearson Chi-squared statistic. Apply your function to the table from (5)(a). You should be able to duplicate the X-squared value (chi-squared) and p-value. Present both.
Shown below are examples of the type of function required. These examples will have to be modified to accomodate the table generated in (5)(a).
chi <- function(x) { # To be used with 2x2 contingency tables that have margins added. # Expected values are calculated. e11 <- x[3,1]x[1,3]/x[3,3] e12 <- x[3,2]x[1,3]/x[3,3] e21 <- x[3,1]x[2,3]/x[3,3] e22 <- x[3,2]x[2,3]/x[3,3] # Value of chi square statistic is calculated. chisqStat <- (x[1,1] - e11)^2/e11 + (x[1,2] - e12)^2/e12 + (x[2,1] - e21)^2/e21 + (x[2,2] - e22)^2/e22 return(list(“chi-squared” = chisqStat, “p-value” = pchisq(chisqStat, 1, lower.tail = F))) }
chisqfun <- function(t) { x <- addmargins(t) e <- matrix(0, nrow = nrow(t), ncol = ncol(t), byrow = T) r <- matrix(0, nrow = nrow(t), ncol = ncol(t), byrow = T) for (i in 1:2) { for (j in 1:2) { e[i,j] = x[nrow(x),j] * x[i,ncol(x)]/x[nrow(x), ncol(x)] r[i,j] = ((x[i,j] - e[i,j])^2)/e[i,j] } } chi <- sum(r) xdf <- nrow(t) - 1 pv <- pchisq(chi, df = xdf, lower.tail = FALSE) return(cat(“Pearson’s Chi-squared test \n”,“Chi sq:”, chi, “; Degree of Freedom :”,xdf," ; P-value :“,pv)) }
chi <- function(x) {
# To be used with 3x2 contingency tables that have margins added.
# Expected values are calculated.
e11 <- x[4,1]*x[1,3]/x[4,3]
e12 <- x[4,2]*x[1,3]/x[4,3]
e21 <- x[4,1]*x[2,3]/x[4,3]
e22 <- x[4,2]*x[2,3]/x[4,3]
e31 <- x[4,1]*x[3,3]/x[4,3]
e32 <- x[4,2]*x[3,3]/x[4,3]
# Value of chi square statistic is calculated.
chisqStat <- (x[1,1] - e11)^2/e11 + (x[1,2] - e12)^2/e12 + (x[2,1] - e21)^2/e21 +
(x[2,2] - e22)^2/e22 + (x[3,1] - e31)^2/e31 + (x[3,2] - e32)^2/e32
return(list("chi-squared" = chisqStat, "p-value" = pchisq(chisqStat, 2, lower.tail = F)))
}
x<-addmargins(contingency)
chi(x)## $`chi-squared`
## [1] 9.086897
##
## $`p-value`
## [1] 0.01063667chisqfun <- function(t) {
x <- addmargins(t)
e <- matrix(0, nrow = nrow(t), ncol = ncol(t), byrow = T)
r <- matrix(0, nrow = nrow(t), ncol = ncol(t), byrow = T)
for (i in 1:3) {
for (j in 1:2) {
e[i,j] = x[nrow(x),j] * x[i,ncol(x)]/x[nrow(x), ncol(x)]
r[i,j] = ((x[i,j] - e[i,j])^2)/e[i,j]
}
}
chi <- sum(r)
xdf <- nrow(t) - 1
pv <- pchisq(chi, df = xdf, lower.tail = FALSE)
return(cat("Pearson's Chi-squared test \\n","Chi sq: ", chi, ";
Degree of Freedom :",xdf," ; P-value :",pv))
}
chisqfun(contingency)## Pearson's Chi-squared test \n Chi sq: 9.086897 ;
## Degree of Freedom : 2 ; P-value : 0.01063667