Instructions

R markdown is a plain-text file format for integrating text and R code, and creating transparent, reproducible and interactive reports. An R markdown file (.Rmd) contains metadata, markdown and R code “chunks”, and can be “knit” into numerous output types. Answer the test questions by adding R code to the fenced code areas below each item. Once completed, you will “knit” and submit the resulting .html file, as well the .Rmd file. The .html will include your R code and the output.

Before proceeding, look to the top of the .Rmd for the (YAML) metadata block, where the title and output are given. Please change title from ‘Programming with R Test #2’ to your name, with the format ‘lastName_firstName.’

If you encounter issues knitting the .html, please send an email via Canvas to your TA.

Each code chunk is delineated by six (6) backticks; three (3) at the start and three (3) at the end. After the opening ticks, arguments are passed to the code chunk and in curly brackets. Please do not add or remove backticks, or modify the arguments or values inside the curly brackets.

Depending on the problem, grading will be based on: 1) the correct result, 2) coding efficiency and 3) graphical presentation features (labeling, colors, size, legibility, etc). I will be looking for well-rendered displays. In the “knit” document, only those results specified in the problem statements should be displayed. For example, do not output - i.e. send to the Console - the contents of vectors or data frames unless requested by the problem. You should be able to display each solution in fewer than ten lines of code.

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 starts from here - There are 5 sections - 75 points total

Section 1: (15 points)
(1) R has probability functions available for use (Kabacoff, Section 5.2.3). Using one distribution to approximate another is not uncommon.

(1)(a) (6 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.

  1. The probability of exactly 0 successes.
dpois(0,lambda = 5)
## [1] 0.006737947
ppois(0,lambda = 5)
## [1] 0.006737947
dbinom(0,100,0.05)
## [1] 0.005920529
pbinom(0,100,0.05)
## [1] 0.005920529
  1. The probability of fewer than 6 successes.
sum(dpois(0:5,lambda = 5))
## [1] 0.6159607
ppois(5,lambda = 5)
## [1] 0.6159607
sum(dbinom(0:5,100,0.05))
## [1] 0.6159991
pbinom(5,100,0.05)
## [1] 0.6159991

(1)(b) (3 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 0.175467370
##  [7] 0.146222808 0.104444863 0.065278039 0.036265577 0.018132789
db<-(dbinom(0:10,100,0.05))
db
##  [1] 0.005920529 0.031160680 0.081181772 0.139575678 0.178142642 0.180017827
##  [7] 0.150014856 0.106025537 0.064870888 0.034901296 0.016715884
par(mfrow=c(1,2))
barplot(dp,main ="Poisson Probability",col="yellow",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.

  1. (3 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.8
varx<-round(((sum((x^2)*prob))-(meanx^2)),digits=2)
varx
## [1] 1.78
  1. (3 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.001
plot(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")

Section 2: (15 points)
(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) (3 points) Load the “faithful” dataset 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.0
str(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)

Wait70<-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      74
erupt3<-subset(Wait70,eruptions<3)
erupt3
##     eruptions waiting
## 211     2.383      71
condprob<-nrow(erupt3)/nrow(Wait70)
condprob
## [1] 0.006060606
  1. (3 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. (1.5 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) (4.5 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   2
colnames(faithful)
## [1] "eruptions" "waiting"
mwait<-matrix(waiting,ncol=2,byrow=TRUE)
plot(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
Section 3: (15 points)
(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) (3 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))
result1<-subset(result,Diet == "1")
boxplot(result1$weight)
result3<-subset(result,Diet == "3")
boxplot(result3$weight)

(3)(b) (3 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.30
Working 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) (3 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) (6 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_post<-mean(post)
mean_weight_gain<-mean_post-mean_pre
mean_weight_gain
## [1] 11.4
n<-length(pre)
df<-n-1
t_lower<-qt(0.025,df)
t_lower
## [1] -2.262157
t_upper<-qt(0.975,df)
t_upper
## [1] 2.262157
std.dev<-sqrt(sum(mean_weight_gain^2)/(n-1))
std.dev
## [1] 3.8
std_error<-std.dev/sqrt(n)
std_error
## [1] 1.201666
lowerci<-mean_weight_gain+(std_error*t_lower)
lowerci
## [1] 8.681644
upperci<-mean_weight_gain+(std_error*t_upper)
upperci
## [1] 14.11836
t.test(pre,post,alternate="two.sided",conf.level = 0.95) #check answers
## 
##  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.3
Section 4: (15 points)
(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) (3 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.3223697
kurtosis1<-kurtosis(Nile)
kurtosis1
## [1] 2.695093
par(mfrow=c(1,2))
qqnorm(Nile,main = "QQ PLOT")
qqline(Nile)
boxplot(Nile,main="Boxplot")

(4)(b) (6 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.7364       42.00156       1764.1312
## sample2 918.5149       20.22883        409.2054

(4)(c) (6 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)

Section 5: (15 points)
(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)(4.5 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] 26
hist(warpbreaks$breaks)
abline(v = median(warpbreaks$breaks),col = "purple",lwd = 2)

number<-ifelse(warpbreaks$breaks<median_breaks,"below","above")
warpbreaks2<-cbind(warpbreaks,number)
summary(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.00
contingency<-table(warpbreaks2$tension,warpbreaks2$number)
contingency
##    
##     above below
##   L    14     4
##   M    10     8
##   H     5    13

(5)(b)(3 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) (7.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 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.01063667
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: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