library(readxl)
library(ggplot2)

Notes on statistical inference made for hadny reference

When sampling is done from a population with mean u and finite standard deviation, the sampling distribution of the sample mean will tend to a normal distribution with mean u and standard deviation sd/sqrt(n) as the sample size n becomes large.

Assumptions:

1. Random Sampling from DGP
2. Stability in DGP
3. If n< 30 => Then assume normal distribution

Congress <- read_excel("Piracy.xlsx", col_names = TRUE)

mean = 30.7833
sd = 1.7862
size = 12
sample_sd = sd/sqrt(size)

lower = mean + qt(0.05, size-1)*sample_sd
higher =mean -qt(0.05, size-1)*sample_sd
lower

## [1] 29.85729

higher

## [1] 31.70931

When sampling is from the same population, using a fixed sample size. The higher the confidence level, wider the confidence band gets. Theoratically at 100% confidence, the band will be infinite wide as that's the only way to be 100% confidenct about the population paramter.

The only way to make confidence interval band narrow is by increasing sample size. ie: get more data!

Assumptions:

1. Random Sampling from DGP
2. Stability in DGP

poll <- read_excel('Poll.xlsx', col_names = TRUE)
table(poll)

## poll
##   1   2 
## 358 407

prop.test can be used for testing the null that the proportions (probabilities of success) in several groups are the same, or that they equal certain given values.

prop.test(358, (358+407), conf.level = .95)

## 
##  1-sample proportions test with continuity correction
## 
## data:  358 out of (358 + 407), null probability 0.5
## X-squared = 3.0118, df = 1, p-value = 0.08266
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.4322160 0.5040576
## sample estimates:
##         p 
## 0.4679739

Performs an exact test of a simple null hypothesis about the probability of success in a Bernoulli experiment. Here hypothesized probability of success is provided as parameter. 0.5 in this case

binom.test(358, (358+407), 0.5, conf.level = .9)

## 
##  Exact binomial test
## 
## data:  358 and (358 + 407)
## number of successes = 358, number of trials = 765, p-value =
## 0.08259
## alternative hypothesis: true probability of success is not equal to 0.5
## 90 percent confidence interval:
##  0.4377759 0.4983554
## sample estimates:
## probability of success 
##              0.4679739

The file WaitTime.xlsx contains data of wait times (in minutes) for a random sample of a bank’s customers. Find a 90% confidence interval for the mean wait time at the bank.

t.test - Performs one and two sample t-tests on vectors of data and also give CI. There are two kinds of hypotheses for a one sample t-test, the null hypothesis and the alternative hypothesis. The alternative hypothesis assumes that some difference exists between the true mean (μ) and the comparison value (m0), whereas the null hypothesis assumes that no difference exists. The comparison value (m0) can be specified with mu parameter below, default is 0. The value of mu does not affect the CI, as CI is the estimation of true mean using leveraging the given sample.

waiting = read_excel('WaitTime.xlsx')
t.test(waiting$WaitTime, mu = 0, conf.level = .90)

## 
##  One Sample t-test
## 
## data:  waiting$WaitTime
## t = 22.057, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 90 percent confidence interval:
##  5.048978 5.871022
## sample estimates:
## mean of x 
##      5.46

t.test(waiting$WaitTime, mu = mean(waiting$WaitTime), conf.level = .90)

## 
##  One Sample t-test
## 
## data:  waiting$WaitTime
## t = 0, df = 99, p-value = 1
## alternative hypothesis: true mean is not equal to 5.46
## 90 percent confidence interval:
##  5.048978 5.871022
## sample estimates:
## mean of x 
##      5.46

An equivalent alternative to t.test() is manual calculation as done below

sample_mean = mean(waiting$WaitTime)
sample_sd = sd(waiting$WaitTime)
sampled_sd = sample_sd/sqrt(length(waiting$WaitTime))

sample_mean+qt(0.05, length(waiting$WaitTime-1))*sampled_sd

## [1] 5.049016

sample_mean-qt(0.05, length(waiting$WaitTime-1))*sampled_sd

## [1] 5.870984

• *Interpretation**: There is a 90% chance that the interval we got captures the true mean waiting time*

Type 1

• Incorrectly rejecting H0 when its actually true
• Hence, you end up accpeting Ha when infact its not true
• Type I error is to falsely infer the existence of something that is not there [conforming to common belief with false information]
• False positive finding
• It is asserting something that is absent, a false hit
• The probability of marking Type I error is called α or Significance-> level
• α is the probability of the study rejecting the null hypothesis, given that it were true
• In statistical hypothesis testing, a result has statistical significance when p < α
• When you have a lower threshold for significance, your probability of concluding there’s an effect when it’s not really there is lower

Type 2

• Failure to reject H0 when its actually false
• Hence, you end up not accepting Ha when infact its true
• Type II error is to falsely infer the absence of something that is present [going against the common belief with false information]
• False negative finding
• It is failing to assert what is present, a miss
• The probability of marking Type II error is called β and and related to the power of a test (which equals 1−β).
• β is the probability of incorrectly concluding no statistical significance
• Statistical power is the likelihood that a study will detect an effect when there is an effect there to be detected.
• Statistical power is the probability that it will reject a false null hypothesis

Example:

• H0 of a fire alarm, is that there is no fire
• Ha of a fire alarm, is that there is a fire
• Type I -> Fire alarm going on indicating a fire when in fact there is no fire
• Type II -> A fire breaking out and the fire alarm does not ring;

The closer is true mean to the hypothetical mean, the higher the value of Beta. We can only control Type 1 error through desired level of significance. Hence, design hypothesis such that the error we want to control is the Type 1 error.

Statisticians Dodge: Since we cannot control Type 2 error. (ie. Accepting Ho when Ha is true). Never draw the conclusion of accepting Ho, instead what we can conclude is that: 'we do not reject Ho) ie: Failed to reject the Null Hypothesis P value is the probability of getting the sample result if Ho were true

The one sample t-test is a statistical procedure used to determine whether a sample of observations could have been generated by a process with a specific mean. Suppose you are interested in determining whether an waiting time at a clinic is less then 5 min. To test this hypothesis, you could collect a sample of waiting times, measure their weights, and compare the sample with a value of five using a one-sample t-test.

The two-sample test is a test performed on the data of two random samples, each independently obtained from a different given population. The purpose of the test is to determine whether the difference between these two populations is statistically significant. One of the model assumptions of the two-sample t-tests for means is that the observations between groups, as well as within groups, are independent. Thus if samples are chosen so that there is some natural pairing, then the two-sample t-test is not appropriate.

Manually calculating T value using the formula $$ t = \frac{(\text{mean}_f - \text{mean}_m) - \text{expected difference}}{SE} \\ ~\\ ~\\ SE = \sqrt{\frac{sd_f^2}{n_f} + \frac{sd_m^2}{n_m}} \\ ~\\ ~\\ \text{where, }~~~df = n_m + n_f - 2 $$ P value of a T statistic can then be calculated as

se = sd(waiting$WaitTime)/sqrt(length(waiting$WaitTime))
t= (mean(waiting$WaitTime)-5)/se
df=length(waiting$WaitTime)-1
pt(t, df)

## [1] 0.9669472

Assumptions:

1. Random Sampling from DGP
2. Stability in DGP
3. If n< 30 => Then assume normal distribution
4. Null hypothesis is true

Test if the waiting time is less than 5 min

Null hypothesis -> true mean is greater then and equal to 5

Alt. hypothesis -> true mean is less than 5

t.test(waiting$WaitTime, alternative = 'less', mu=5)

## 
##  One Sample t-test
## 
## data:  waiting$WaitTime
## t = 1.8582, df = 99, p-value = 0.9669
## alternative hypothesis: true mean is less than 5
## 95 percent confidence interval:
##      -Inf 5.871022
## sample estimates:
## mean of x 
##      5.46

The likely hood of getting sample statistic as extreme as 1.85 is 0.966 if the true mean waiting time is greater then or equal to 5

The P Value is the probability of obtaining a result as extreme as, or more extreme than, the result actually obtained when the null hypothesis is true

The significance level and confidence level are the complementary portions in the normal distribution.

Test normality for generated numbers

Assumptions:

1. Random Sampling from DGP
2. Stability in DGP

seq <- rnorm(50,5,1)
hist(seq)

qqnorm(seq)
qqline(seq, col = "blue")

Null hypothesis -> distribution is normal

Alt. hypothesis -> distribution is not normal

shapiro.test(seq)

## 
##  Shapiro-Wilk normality test
## 
## data:  seq
## W = 0.96493, p-value = 0.1429

P-value being so high, we can conclude that: We fail to reject null hypothesis. Hence, our assumption of distribution being normal is not rejected.

Definition of P-Value: Given a null hypothesis and sample evidence of size n, the p-value is the probability of getting a sample evidence that is equally or more unfavorable to the null hypothesis while the null hypothesis is actually true.

The Chi Square distribution is the distribution of the sum of squared standard normal deviates. The degrees of freedom of the distribution is equal to the number of standard normal deviates being summed

Probability of getting X-square value of less then or equal to 5 is

pchisq(5,1)

## [1] 0.9746527

X-square value associated with cumulative probability of 0.9746527 is

qchisq(0.9746527,1)

## [1] 5.000001

The chi-squared test is used to determine whether there is a significant difference between the expected frequencies and the observed frequencies in one or more categories.

Chi-Square goodness of fit test is a non-parametric test that is used to find out how the observed value of a given phenomena is significantly different from the expected value. In Chi-Square goodness of fit test, the term goodness of fit is used to compare the observed sample distribution with the expected probability distribution.

Null hypothesis -> There is no significant difference between the observed and the expected value.

Alt. hypothesis -> There is a significant difference between at least one of them

Assumptions:

1. Random Sampling from DGP
2. Stability in DGP
3. Expected value for each group is at least 5

gas<- read_excel("Gasoline.xlsx", col_names=TRUE)
head(gas)

## # A tibble: 6 x 3
##   Owner `Second-last`  Last
##   <dbl>         <dbl> <dbl>
## 1     1             3     3
## 2     2             4     3
## 3     3             2     2
## 4     4             3     3
## 5     5             4     3
## 6     6             2     4

table(gas$`Second-last`, gas$Last)

##    
##      1  2  3  4
##   1 47 14 27 14
##   2 18 63 33 13
##   3 26 48 59 16
##   4 30 23 29 53

result <- chisq.test(table(gas$`Second-last`, gas$Last))
result

## 
##  Pearson's Chi-squared test
## 
## data:  table(gas$`Second-last`, gas$Last)
## X-squared = 114.12, df = 9, p-value < 2.2e-16

Check expected values are greater then equal to 5

result$expected

##    
##            1        2        3        4
##   1 24.05848 29.42690 29.42690 19.08772
##   2 29.95517 36.63938 36.63938 23.76608
##   3 35.14425 42.98635 42.98635 27.88304
##   4 31.84211 38.94737 38.94737 25.26316

Cell-by-cell contributions to observed Chi-Square value

(result$residuals)^2

##    
##              1          2          3          4
##   1 21.8764183  8.0874729  0.2001518  1.3561017
##   2  4.7713302 18.9654562  0.3614979  4.8770563
##   3  2.3792598  0.5847585  5.9655403  5.0642490
##   4  0.1065681  6.5298009  2.5406117 30.4527412

Null hypothesis -> No Relationship

Alt. hypothesis -> Relationship

Looking at relationships between variables to understand and and predict what is going on.

Ho:No relationship exists; probabilistic independence

Ha: Relationship exists; probabilistic dependence

Ex. Transit Railroads is interested in the relationship between travel distance and the ticket class purchased. A random sample of 200 passengers is taken.

transit<-read_excel("Transit.xlsx", na="NA", col_names = TRUE)
table(transit)

##         Distance
## Class    1-100 101-200 201-300 301-400 401-500
##   First      6       8      15      21      10
##   Second    14      16      17      14       6
##   Third     21      18      16      12       6

result <- chisq.test(table(transit))
result

## 
##  Pearson's Chi-squared test
## 
## data:  table(transit)
## X-squared = 15.923, df = 8, p-value = 0.0435

# cell-by-cell contributions to observed Chi-Square value: ((O-E)^2)/E
result$residuals

##         Distance
## Class         1-100    101-200    201-300    301-400    401-500
##   First  -1.7963377 -1.2959032  0.1581139  1.8375516  1.3234482
##   Second  0.0715042  0.5145295  0.2294271 -0.4397685 -0.5046460
##   Third   1.5600514  0.6819306 -0.3631420 -1.2446101 -0.7163714

F - Distribution: Generally arises from a statistic that involves a ratio of variances. pf & qf functions in R for handling F-distribution

F statistic is the value we receive when we run an ANOVA test on different groups to understand the differences between them. The F statistic is given by the ratio of between group variability to within group variability

If there are two predictors, then it becomes Twoway ANOVA Oneway ANOVA is a test of relationship between a interval level response variable and nominal level predictor variable. Or can also be expressed as a test of multiple means.

Assumptions:

1. Random Sampling from DGP
2. Stability in DGP
3. For each level of predictor variable, if n < 30 then we assume normal distrubtion
4. Standard deviation is same for all the groups (Evaluate with Leven's test)

Ho:No relationship exists; Distribution is identical across all categories

Ha: Relationship exists; Difference in means across at-least two categories

chol<- read_excel("Cholesterol.xlsx", col_names=TRUE)
head(chol)

## # A tibble: 6 x 2
##   Drug  CholReduction
##   <chr>         <dbl>
## 1 A                22
## 2 A                31
## 3 A                19
## 4 A                27
## 5 A                25
## 6 A                18

fit <-  aov(CholReduction ~ Drug, data=chol)
summary(fit)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Drug         2 2152.1  1076.1   40.79 8.59e-07 ***
## Residuals   15  395.7    26.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

DF = (k-1), (n-k) = 2, 15 Mean Sq, gives MST (above) and MSE (below) Sum Sq, gives SST (above) and SSE (below) The probability of getting a sample statistic as or more extreme then 40.79 is the P value. Which is extremely low. Hence, there is a significant relationship between means of at least 2 drugs.

boxplot(CholReduction ~ Drug, data=chol)

Use Tukey HSD test to determine which drugs have relationship.

Studentised range distribution is used to test each pairwise difference

Assumptions: Same as ANOVA

1. Random Sampling from DGP
2. Stability in DGP
3. For each level of predictor variable, in n < 30 then we assume normal distribution
4. Standard deviation is same for all the groups (Evaluate with Leven's test)

Ho:No relationship exists: μi = μj for each i, j being checked

Ha: Relationship exists: μi ≠ μj

TukeyHSD(fit, conf.level = .9)

##   Tukey multiple comparisons of means
##     90% family-wise confidence level
## 
## Fit: aov(formula = CholReduction ~ Drug, data = chol)
## 
## $Drug
##          diff        lwr       upr     p adj
## B-A  15.50000   8.916894  22.08311 0.0002841
## C-A -11.16667 -17.749773  -4.58356 0.0049961
## C-B -26.66667 -33.249773 -20.08356 0.0000006

Ex. here. B is greater then A buy 15.4 There is a probability of 0.00028 or more to get a sample statistic of 15.5 if Ho where true Hence, all three of them differ significantly

Whenever we wish to determine if variances are all equal or not across different groups

Assumptions:

1. Random Sampling from DGP
2. Stability in DGP

Ho: σ1 = σ2 = … = σk ; k independent samples/sets of responses

Ha: At least one pair is not equal

library('car')

## Loading required package: carData

leveneTest(chol$CholReduction, chol$Drug)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.0882 0.9161
##       15

Assumptions:

1. Random Sampling from DGP
2. Stability in DGP
3. εi ~ Normal (Mean 0, SD σ) constant, over all possible values of the predictors. We can use the following to evaluate this assumption
   • Standardized Residuals are normally distributed
   • Mean 0, SD fixed across predictors. Plot St. residual plots for individual predictors (Xi) and predicted values (Y hat)

We are interested in predicting the amount of carbohydrates (in grams) for a menu item based on its calorie content (measured in 100s).

starbucks <- read_excel('starbucks.xlsx')
head(starbucks)

## # A tibble: 6 x 7
##   item                        calories   fat  carb fiber protein type  
##   <chr>                          <dbl> <dbl> <dbl> <dbl>   <dbl> <chr> 
## 1 8-Grain Roll                     3.5     8    67     5      10 bakery
## 2 Apple Bran Muffin                3.5     9    64     7       6 bakery
## 3 Apple Fritter                    4.2    20    59     0       5 bakery
## 4 Banana Nut Loaf                  4.9    19    75     4       7 bakery
## 5 Birthday Cake Mini Doughnut      1.3     6    17     0       0 bakery
## 6 Blueberry Oat Bar                3.7    14    47     5       6 bakery

attach(starbucks)
linefit <- lm(carb ~ calories, data = starbucks)
plot(starbucks$calories, starbucks$carb)
abline(linefit,lty=2, col="blue")

Ho:No relationship exists; β1 = β2 =…=βk = 0

Ha: Relationship exists; at least one predictor βi ≠ 0

ANOVA is used to test this

anova(linefit)

## Analysis of Variance Table
## 
## Response: carb
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## calories   1  9486.4  9486.4  62.772 1.673e-11 ***
## Residuals 75 11334.3   151.1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ho:No relationship exists; βi = 0, 1 ≤ i ≤ k

Ha: Relationship exists; βi ≠ 0

For each predictor two tailed T - test is used to test this

summary(linefit)

## 
## Call:
## lm(formula = carb ~ calories, data = starbucks)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.477  -7.476  -1.029  10.127  28.644 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    8.944      4.746   1.884   0.0634 .  
## calories      10.603      1.338   7.923 1.67e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.29 on 75 degrees of freedom
## Multiple R-squared:  0.4556, Adjusted R-squared:  0.4484 
## F-statistic: 62.77 on 1 and 75 DF,  p-value: 1.673e-11

Standardized residuals plot

linefit.stres <- rstandard(linefit)

plot(starbucks$calories, linefit.stres, pch = 16, main = "Standardized Residual Plot", xlab = "Calories", ylab = "Standardized Residuals
")
abline(0,0, lty=2, col="blue")

Normal probability plot (QQPlot)

qqnorm(linefit.stres, main = "Normal Probability Plot", xlab = "Normal Scores", ylab = "Standardized Residuals")
qqline(linefit.stres, col = "blue")

Shapiro Wilk Normality test

shapiro.test(linefit.stres)

## 
##  Shapiro-Wilk normality test
## 
## data:  linefit.stres
## W = 0.99032, p-value = 0.832

plot(linefit$fitted.values, linefit.stres, pch = 16, main = "Standardized Residual Plot", xlab = "Fitted Carbs", ylab = "Standardized Residuals")
abline(0,0, lty=2, col="red")

plot(starbucks$calories, linefit.stres, pch = 16, main = "Standardized Residual Plot", xlab = "Fitted Chol", ylab = "Standardized Residuals")
abline(0,0, lty=2, col="red")

confint(linefit, level=.9)

##                  5 %     95 %
## (Intercept) 1.039446 16.84767
## calories    8.374277 12.83190

predict(linefit,data.frame(calories=4.5), interval ="confidence", level = .80)

##        fit      lwr      upr
## 1 56.65746 54.01528 59.29964

predict(linefit, data.frame(calories=4.5), interval="predict", level = .80)

##        fit     lwr      upr
## 1 56.65746 40.5449 72.77002