3. Probability Analysis Part 1.Rmd

---
title: 'Gowani_Ali'
output:
  html_document: default
---

```{r setup, include=FALSE}
# DO NOT ADD OR REVISE CODE HERE
knitr::opts_chunk$set(echo = TRUE, eval = TRUE)
```


-----

### 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.

(i) The probability of exactly 0 successes.

```{r test1ai}

dpois(0, lambda = 5)

ppois(0, lambda = 5)

dbinom(0, 100, 0.05)

pbinom(0, 100, 0.05)

```

(ii) The probability of fewer than 6 successes.

```{r test1aii}

sum(dpois(0:5,lambda = 5))

ppois(5, lambda = 5)

sum(dbinom(0:5, 100, 0.05))

pbinom(5, 100, 0.05)

```

(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.

```{r test1b}

library(ggplot2)
library(gridExtra)

plot_dpois <- dpois(0:10, lambda = 5)
plot_dbinom <- dbinom(0:10, 100, 0.05)

ggplot_dpois <- ggplot(data = NULL, aes(x = 0:10, y = plot_dpois)) + 
                     geom_col(fill = "#56B4E9", alpha = .6) +
                     labs(subtitle="Poisson Probabilities", 
                         y = "Poisson Probabilities", 
                         x = "Range (0 - 10)", 
                         title = "R Assignment #2, Section 1.B", 
                         caption = "2019SU_MSDS_401-DL_SEC55 - Prof. Srinivasan") +
                     scale_x_discrete(name = "Range (0 - 10)", limits = 0:10)
  
ggplot_pbinom <- ggplot(data = NULL, aes(x = 0:10, y = plot_dbinom)) + 
                     geom_col(fill = "#E69F00", alpha = .6) +
                     labs(subtitle="Binomial Probabilities", 
                        y = "Binomial Probabilities", 
                        x = "Range (0 - 10)", 
                        title = "R Assignment #2, Section 1.B", 
                        caption = "2019SU_MSDS_401-DL_SEC55 - Prof. Srinivasan") +
                     scale_x_discrete(name = "Range (0 - 10)", limits = 0:10)

grid.arrange(ggplot_dpois, ggplot_pbinom, ncol = 2, nrow = 1)

```

(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.  

(i) (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. 

```{r test1ci}

dataframe_ev <- data.frame("Outcome" = c(0:6), "Probabilities" = c(0.215, 0.230, 0.240, 0.182, 0.130, 0.003, 0.001))

dataframe_ev

ev <- sum(dataframe_ev$Outcome * dataframe_ev$Probabilities)

print(round(ev, digits = 2))

print(round(var(dataframe_ev$Probabilities), digits = 2))


```

(ii) (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.

```{r test1cii}

cumsum_probs <- cumsum(dataframe_ev$Probabilities)

cumsum_probs_median <- median(cumsum_probs)
   
summary(cumsum_probs)

cat("Median: ", cumsum_probs_median)

ggplot(data = dataframe_ev, aes(dataframe_ev$Outcome)) + 
  geom_point(aes(y = cumsum_probs, color = "Cummulative Probability")) + 
  labs(subtitle="Cummulative Probabilities",  
      y = "Cummulative Probabilities", 
      x = "Range (0 - 10)", 
      title = "R Assignment #2, Section 1.C.2", 
      caption = "2019SU_MSDS_401-DL_SEC55 - Prof. Srinivasan") +
   scale_x_discrete(name = "Outcomes (0 - 6)", limits = 0:6) +
   geom_abline(slope = 0, intercept = cumsum_probs_median, color = "#E69F00") +
   scale_fill_discrete(guide=FALSE) +
   geom_text(size=3, aes(0,.9), label = "Median")

```

##### 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.

```{r test2a}

library(mosaic)

data(faithful)

summary(faithful)

plot_faithful <- ggplot(faithful, aes(x = faithful$waiting)) +
                     geom_histogram(fill = "red", alpha = .6, bins=12) +
                     labs(y = "Frequency", 
                          x = "Waiting Times Between Eruptions (in minutes)", 
                          subtitle = "R Assignment #2, Section 2.A")

plot_faithful

print(tally(~faithful$eruptions < 3 & faithful$waiting > 70, data = faithful, margins=TRUE), digits = 2)

point_370 <- round(count(faithful$eruptions < 3 & faithful$waiting > 70) / count(faithful$waiting > 70) * 100, digits = 2)

probability_faithful <- (count(faithful$eruptions < 3 & faithful$waiting > 70) / count(faithful))

cat("Probability of eruption less than 3.0 minutes, if the waiting time exceeds 70 minutes: ", paste(round(count(faithful$eruptions < 3 & faithful$waiting > 70) / count(faithful$waiting > 70) * 100, digits = 2), "%", sep = ""))

```

(i) (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. 

```{r test2ai}

data(faithful)

erupt3_wait70 <- which(faithful$eruptions < 3.0 & faithful$waiting > 70)

print(faithful[erupt3_wait70,])

ggplot(data = faithful, aes(x = faithful$waiting)) + 
   geom_point(aes(x = faithful$waiting[erupt3_wait70], faithful$eruptions[erupt3_wait70], alpha = 1/10), size = 4, alpha = 1/10, col="orange") +
  geom_point(aes(y = faithful$eruptions, color = "Eruption")) + 
  labs(subtitle="R Assignment", 
       y = "eruption minutes", 
       x = "wait times (minutes)",
       title = "R Assignment #2, Section 2.A.I", 
       caption = "2019SU_MSDS_401-DL_SEC55 - Prof. Srinivasan") +
   geom_hline(aes(yintercept=3), color = "orange", alpha = .6) +
   geom_vline(aes(xintercept=70), color = "orange", alpha = .6) + 
   scale_color_discrete(name = "Event:") +
   geom_text(size=3, aes(80,2.2), label = "Eruption where wait time \n exceeds 70 minutes and duration \n is less than 3.0 minutes.", colour = "orange")

```

(ii) (1.5 point) What does the plot suggest about the relationship between eruption time and waiting time?

***Answer: In general, there seems to be a positive correlation between the duration of the eruption and the wait times associated with the eruption. For example, the longer a person waits for the eruption to occur, then there is a greater likleyhood that the person will see a longer eruption. However, looking at the chart and in particular the horizontal and verticle lines, we can strongly suggest that if a person waits 70 minutes for an eruption then that eruption should last for more than 3.0 minutes. Out of 272 observations in our dataset, there was 1 example where this did not hold true.  ***

-----

(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.

```{r test2b}

new_faithful <- data.frame(faithful)

new_faithful_odd <- new_faithful[c(TRUE, FALSE),]     # for odd rows
new_faithful_even <- new_faithful[c(FALSE, TRUE),]     # for even rows

odd_even_faithful <- cbind(first=new_faithful_odd$waiting, second=new_faithful_even$waiting)
odd_even_faithful <- data.frame(odd_even_faithful)

ggplot(data = odd_even_faithful, aes(x = odd_even_faithful$first)) + 
  geom_point(aes(y = odd_even_faithful$second, color = "Eruption")) + 
  labs(subtitle="Are wait times between eruptions independent?", 
       y = "second eruption", 
       x = "first eruption",
       title = "R Assignment #2, Section 2.B", 
       caption = "2019SU_MSDS_401-DL_SEC55 - Prof. Srinivasan")

```

(2)(c) (2) Test the hypothesis of independence with a two-sided test at the 5% level using the Kendall correlation coefficient.  

```{r test2c}

kendall_test <- cor.test(odd_even_faithful$first, odd_even_faithful$second,  method="kendall", Cor.coeff = .05, alternative = "two.sided")

kendall_test

```

##### 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.

```{r test3}

# 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) 

```

##### 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". 

```{r test3a}

boxplot_Diet1 <- ggplot(data = subset(result, Diet == "1"), aes(y = weight, x = Diet)) +
                     stat_boxplot(geom ='errorbar', linetype = 1, width = 0.25) +
                     geom_boxplot(fill = "#56B4E9", alpha = .9) + 
                     labs(y = "Weight", 
                     x = "Diet 1", 
                     subtitle = "Boxplot of weight where Diet is 1")


boxplot_Diet3 <- ggplot(data = subset(result, Diet == "3"), aes(y = weight, x = Diet)) +
                     stat_boxplot(geom ='errorbar', linetype = 1, width = 0.25) +
                     geom_boxplot(fill = "#E69F00", alpha = .9) + 
                     labs(y = "Weight", 
                     x = "Diet 3", 
                     subtitle = "Boxplot of weight where Diet is 3")

grid.arrange(boxplot_Diet1, boxplot_Diet3, ncol = 2, nrow = 1)
 
```

(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().

```{r test3b}

weight_Diet1 <- result$weight[result$Diet == "1"]
weight_Diet3 <- result$weight[result$Diet == "3"]

t.test(x = weight_Diet1, y = weight_Diet3, alternative = "two.sided", mu = 0, 
       paired = FALSE, var.equal = FALSE, conf.level = 0.95)

```

##### 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.

```{r test3paired}

# 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)

```

(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.  

```{r test3c}

prepost_data <- cbind(pre, post)
df.prepost_data <- data.frame(prepost_data)

ggplot(data = df.prepost_data, aes(df.prepost_data$pre)) + 
  geom_point(aes(y = df.prepost_data$post, color = df.prepost_data$post)) + 
  ggtitle("Abalones: WHOLE versus SHUCK") +
  labs(subtitle="Variable Post Weight as function of Pre Weight.", 
       y = "Post Weight", 
       x = "Pre Weight", 
       title = "R Assignment #2, Section 3.C", 
       caption = "2019SU_MSDS_401-DL_SEC55 - Prof. Srinivasan") + 
  geom_abline(slope = 1, intercept = 0, col = "firebrick", alpha = .6) +
  scale_color_gradient(low = 'darkgreen', high = 'gold3', name = "Weight")

```

(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.

```{r test3d}

pre_mean <- mean(df.prepost_data$pre)
post_mean <- mean(df.prepost_data$post)

predataset <- df.prepost_data$pre
postdataset <- df.prepost_data$post

var_predataset <- var((predataset)/(length(predataset)))
var_postdataset <- var((postdataset)/(length(postdataset)))

prepost_mean <- (sum(df.prepost_data$post) - sum(df.prepost_data$pre)) / length(df.prepost_data$pre)
prepost_mean
prepost_n <- length(df.prepost_data$pre)
prepost_n
prepost_df <- prepost_n - 1
prepost_df

lower_t <- qt(0.025, prepost_df)
upper_t <- qt(0.975, prepost_df)

prepost_sd <- sqrt(sum(prepost_mean^2)/(prepost_n-1))
prepost_sd
prepost_sd_error <- prepost_sd / sqrt(prepost_n)
prepost_sd_error

prepost_confint_lower <- prepost_mean + (prepost_sd_error * lower_t)
prepost_confint_upper <- prepost_mean + (prepost_sd_error * upper_t)

p.value = dt(upper_t, df=prepost_df)
p.value

t.value <- (mean(df.prepost_data$post) - mean(df.prepost_data$pre)) / sqrt((var((df.prepost_data$post) / prepost_n)) + var((df.prepost_data$pre) / prepost_n))
t.value

pt.value <- pt(t.value,df=prepost_df)
pt.value

```

##### 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.

```{r test4}

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.

```{r test4a}

library(moments) 

nile_skewness <- skewness(Nile)
nile_skewness

nile_kurtosis <- kurtosis(Nile)
nile_kurtosis

df.nile <- data.frame(Nile)

nile_qq <- ggplot(df.nile, aes(sample = Nile, alpha = .6)) +
             stat_qq(show.legend = FALSE, color = "#00AFBB", fill = "#C4961A") + 
             stat_qq_line(show.legend = FALSE, color = "#C4961A", fill = "#C4961A") +
             labs(y = "Flow", 
                  x = "Theoretical Quantiles", 
                  subtitle = "QQ Plot & QQ Line: Flow Data")

nile_box <- ggplot(df.nile, aes(y = Nile)) +
               stat_boxplot(geom ='errorbar', linetype = 1, width = 0.25, color = "#C4961A") +
               geom_boxplot(alpha = .7, fill = "#C4961A", color = "#C4961A") +
               labs(y = "Flow", 
                    x = "Nile",
                    subtitle = "Box Plot: Flow Data") + 
               theme(axis.text.x=element_blank())

grid.arrange(nile_qq, nile_box, ncol = 2, nrow = 1)

```

(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.

```{r test4b}

library(kableExtra)

set.seed(124)
sample1 <- rep(1:1000, 0)

for (i in 1:1000) {
   sample1[i] <- mean(sample(Nile, n = 16, replace = TRUE))
}

set.seed(127)
sample2 <- rep(1:1000, 0)

for (i in 1:1000) {
  sample2[i] <- mean(sample(Nile, n = 64, replace = TRUE))
}

mean_sample1 <- mean(sample1)
mean_sample2 <- mean(sample2)

sd_sample1 <- sd(sample1)
sd_sample2 <- sd(sample2)

var_sample1 <- var(sample1)
var_sample2 <- var(sample2)

trib_nile <-  tribble(~"Sample Type",    ~Mean,   ~"Sample SD", ~"Sample Variance",
                        "Sample 1", mean_sample1, sd_sample1, var_sample1,
                        "Sample 1", mean_sample1, sd_sample2, var_sample1,
                        )
trib_nile

kable(trib_nile, format = "html", caption = "Nile Sample Statistics") %>%
      kable_styling(bootstrap_options = "striped",
                    full_width = F,
                    font_size = 12,
                    position = "left") %>%
      add_footnote(c("R Assignment #2, Section 4.B"))

```

(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.** 

```{r test4c}

par(mfrow=c(1,2))

hist(sample1, 
     main="Sample 1: Histogram with Curve", 
     xlab="Sample 1",
     freq = FALSE, 
     xlim = c(750, 1040), 
     ylim = c(0, 0.025))
curve(dnorm(x, 
            mean = mean(sample1), 
            sd = sd(sample1)), 
      add=TRUE, 
      col="blue")

hist(sample2, 
     main="Sample 2: Histogram with Curve", 
     xlab="Sample 2",
     freq = FALSE, 
     xlim = c(750, 1040), 
     ylim = c(0, 0.025))
curve(dnorm(x, 
            mean = mean(sample2), 
            sd = sd(sample2)), 
      add=TRUE, 
      col="blue")

df_sample1 <- data.frame(sample1)
df_sample2 <- data.frame(sample2)

plot_sample1 <- ggplot(data = df_sample1, aes(x = df_sample1$sample1, y = ..density..,)) +
                  geom_histogram(fill = "red", alpha = .6, bins=11) +
                  labs(y = "Density", 
                       x = "Sample 1", 
                       subtitle = "Sample 1: Histogram with Curve")

plot_sample2 <- ggplot(data = df_sample2, aes(x = df_sample2$sample2, y = ..density..,)) +
                  geom_histogram(fill = "red", alpha = .6, bins=11) +
                  labs(y = "Density", 
                       x = "Sample 2", 
                       subtitle = "Sample 2: Histogram with Curve")

grid.arrange(plot_sample1, plot_sample2, ncol = 2, nrow = 1)

```

-----

##### 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.

```{r test5a}

data(warpbreaks)

str(warpbreaks)

median_breaks <- median(warpbreaks$breaks)

ggplot(warpbreaks, aes(x = warpbreaks$breaks)) +
      geom_histogram(fill = "#E69F00", alpha = .6, bins = 5) +
      labs(y = "Frequency", 
           x = "Breaks", 
           subtitle = "Warp Breaks") +
      geom_vline(xintercept = median_breaks, color = "#D55E00", size = 1.5) +
      scale_fill_discrete(guide=FALSE) +
      geom_text(size=3, aes(29,.5), label = "Median")

number <- ifelse(
            warpbreaks$breaks<median_breaks,
            "Below",
            "Above")

warpbreaks <- cbind(warpbreaks, number)

summary(warpbreaks)

contingency_table <- table(warpbreaks$tension, warpbreaks$number)
contingency_table

contingency_table_margins <- addmargins(table(warpbreaks$tension, warpbreaks$number))

kable(contingency_table_margins, format = "html", caption = "Contingency Table") %>%
      kable_styling(bootstrap_options = "striped",
                    full_width = F,
                    font_size = 12,
                    position = "left") %>%
      add_footnote(c("R Assignment #2, Section 5.A"))

```

(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. 

```{r test5b}

chisq.test(contingency_table, correct = FALSE)

print("Based on the chisq.test, we get the chi-squared value of 9.0869 and we get a p-value of 0.01064. The p-value (0.01), is less than the significance level of 0.05. Therefore, we can reject the null hypothesis as the two variables (tensions and breaks) are dependent.")

```


(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).  

```{r test5c}

chi <- function(x) {
   # To be used with 2x2 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)))
}

chi_x <- addmargins(contingency_table)
chi(chi_x)


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_table)

```