Test of Hypotheses

Author

Marcio Diniz

Published

November 7, 2022

Libraries

# install.packages("BSDA")
library(BSDA)

Loading required package: lattice


Attaching package: 'BSDA'

The following object is masked from 'package:datasets':

    Orange

library(janitor)


Attaching package: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test

library(tidyverse)

── Attaching packages
───────────────────────────────────────
tidyverse 1.3.2 ──

✔ ggplot2 3.3.6     ✔ purrr   0.3.4
✔ tibble  3.1.8     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.1
✔ readr   2.1.2     ✔ forcats 0.5.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

#install.packages("nortest")
library(nortest)
#install.packages("lawstat")
library(lawstat)
#install.packages("nparcomp")
library(nparcomp)

Loading required package: multcomp
Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data
Loading required package: MASS

Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select


Attaching package: 'TH.data'

The following object is masked from 'package:MASS':

    geyser

Datasets

dataset <- read_csv(file = "data/troponin.csv") %>% 
  clean_names(case = "old_janitor")

Rows: 369 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Gender
dbl (2): ID, TnI Conc

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ggplot(dataset, aes(x = tni_conc, y = stat(density))) +
  geom_histogram() + 
  labs(x = "Troponin I", y = "Density")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

dt <- dataset %>%
  filter(tni_conc < 900)

One sample histogram

ggplot(dt, aes(x = tni_conc, y = after_stat(density))) +
  geom_histogram() + 
  labs(x = "Troponin I", y = "Density")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

One Sample z-test

# Testing using the R-package BSDA
z.test(x = dt$tni_conc, alternative = "two.sided", mu = 2, sigma.x = 1)


    One-sample z-Test

data:  dt$tni_conc
z = -11.265, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 2
95 percent confidence interval:
 1.310603 1.514943
sample estimates:
mean of x 
 1.412773

z.test(x = dt$tni_conc, alternative = "greater", mu = 2, sigma.x = 1)


    One-sample z-Test

data:  dt$tni_conc
z = -11.265, p-value = 1
alternative hypothesis: true mean is greater than 2
95 percent confidence interval:
 1.327029       NA
sample estimates:
mean of x 
 1.412773

z.test(x = dt$tni_conc, alternative = "less", mu = 2, sigma.x = 1)


    One-sample z-Test

data:  dt$tni_conc
z = -11.265, p-value < 2.2e-16
alternative hypothesis: true mean is less than 2
95 percent confidence interval:
       NA 1.498517
sample estimates:
mean of x 
 1.412773

# Calculating the Z-score
z_score <- (mean(dt$tni_conc) - 2)/sqrt(1/nrow(dt))

# Calculating p value
pnorm(z_score, mean = 0, sd = 1)

[1] 9.771756e-30

p_value <- 2*(1 - pnorm(z_score, mean = 0, sd = 1))

# Variance of our data
sd(dt$tni_conc)

[1] 3.571236

sd(dt$tni_conc)^2

[1] 12.75372

var(dt$tni_conc)

[1] 12.75372

One Sample t-test

### t test

# Testing
t.test(x = dt$tni_conc, alternative = "two.sided", mu = 2)


    One Sample t-test

data:  dt$tni_conc
t = -3.1544, df = 367, p-value = 0.001741
alternative hypothesis: true mean is not equal to 2
95 percent confidence interval:
 1.046692 1.778854
sample estimates:
mean of x 
 1.412773

t.test(x = dt$tni_conc, alternative = "greater", mu = 2)


    One Sample t-test

data:  dt$tni_conc
t = -3.1544, df = 367, p-value = 0.9991
alternative hypothesis: true mean is greater than 2
95 percent confidence interval:
 1.105787      Inf
sample estimates:
mean of x 
 1.412773

t.test(x = dt$tni_conc, alternative = "less", mu = 2)


    One Sample t-test

data:  dt$tni_conc
t = -3.1544, df = 367, p-value = 0.0008706
alternative hypothesis: true mean is less than 2
95 percent confidence interval:
    -Inf 1.71976
sample estimates:
mean of x 
 1.412773

# Calculating the T-score
t_score <- (mean(dt$tni_conc) - 2)/sqrt(var(dt$tni_conc)/nrow(dt))

# p-value for the two-sided test
pt(t_score, df = (nrow(dt) - 1))

[1] 0.0008705882

p_value <- 2*pt(t_score, df = (nrow(dt) - 1))

Two sample histogram

ggplot(dt, aes(x = tni_conc, y = stat(density), fill = gender)) +
  geom_histogram() + 
  facet_grid(. ~ gender, labeller = labeller(gender = label_both)) + 
  theme(legend.position = "none") +
  labs(x = "Troponin I", y = "Density")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Student t-test for two samples with equal variance

# t test with var.equal = TRUE
t.test(tni_conc ~ gender, 
       alternative = "two.sided", var.equal = TRUE,
       paired = FALSE, data = dt)


    Two Sample t-test

data:  tni_conc by gender
t = -2.5462, df = 366, p-value = 0.0113
alternative hypothesis: true difference in means between group F and group M is not equal to 0
95 percent confidence interval:
 -1.6678834 -0.2142782
sample estimates:
mean in group F mean in group M 
      0.9371183       1.8781991

t.test(tni_conc ~ gender, 
       alternative = "greater", var.equal = TRUE,
       paired = FALSE, data = dt)


    Two Sample t-test

data:  tni_conc by gender
t = -2.5462, df = 366, p-value = 0.9944
alternative hypothesis: true difference in means between group F and group M is greater than 0
95 percent confidence interval:
 -1.550558       Inf
sample estimates:
mean in group F mean in group M 
      0.9371183       1.8781991

t.test(tni_conc ~ gender, 
       alternative = "less", var.equal = TRUE,
       paired = FALSE, data = dt)


    Two Sample t-test

data:  tni_conc by gender
t = -2.5462, df = 366, p-value = 0.005649
alternative hypothesis: true difference in means between group F and group M is less than 0
95 percent confidence interval:
       -Inf -0.3316031
sample estimates:
mean in group F mean in group M 
      0.9371183       1.8781991

dt %>% group_by(gender) %>% summarize(var = var(tni_conc))

# A tibble: 2 × 2
  gender   var
  <chr>  <dbl>
1 F       2.38
2 M      22.5

Welch t-test for two samples with unequal variance

t.test(tni_conc ~ gender, 
       alternative = "two.sided", var.equal = FALSE,
       paired = FALSE, data = dt)


    Welch Two Sample t-test

data:  tni_conc by gender
t = -2.5688, df = 224.42, p-value = 0.01085
alternative hypothesis: true difference in means between group F and group M is not equal to 0
95 percent confidence interval:
 -1.6630143 -0.2191473
sample estimates:
mean in group F mean in group M 
      0.9371183       1.8781991

t.test(tni_conc ~ gender, 
       alternative = "greater", var.equal = FALSE,
       paired = FALSE, data = dt)


    Welch Two Sample t-test

data:  tni_conc by gender
t = -2.5688, df = 224.42, p-value = 0.9946
alternative hypothesis: true difference in means between group F and group M is greater than 0
95 percent confidence interval:
 -1.546177       Inf
sample estimates:
mean in group F mean in group M 
      0.9371183       1.8781991

t.test(tni_conc ~ gender, 
       alternative = "less", var.equal = FALSE,
       paired = FALSE, data = dt)


    Welch Two Sample t-test

data:  tni_conc by gender
t = -2.5688, df = 224.42, p-value = 0.005427
alternative hypothesis: true difference in means between group F and group M is less than 0
95 percent confidence interval:
       -Inf -0.3359845
sample estimates:
mean in group F mean in group M 
      0.9371183       1.8781991

Normality tests

ggplot(dt, aes(sample = (tni_conc - mean(tni_conc))/sd(tni_conc))) + stat_qq() + stat_qq_line()

ggplot(dt, aes(sample = (tni_conc - mean(tni_conc))/sd(tni_conc))) + stat_qq() + stat_qq_line() +
  facet_grid(. ~ gender)

# Shapiro-Francia test
sf.test(dt$tni_conc[dt$gender == "F"])


    Shapiro-Francia normality test

data:  dt$tni_conc[dt$gender == "F"]
W = 0.48815, p-value < 2.2e-16

sf.test(dt$tni_conc[dt$gender == "M"])


    Shapiro-Francia normality test

data:  dt$tni_conc[dt$gender == "M"]
W = 0.26448, p-value < 2.2e-16

# Anderson-Darling test

ad.test(dt$tni_conc[dt$gender == "F"])


    Anderson-Darling normality test

data:  dt$tni_conc[dt$gender == "F"]
A = 30.439, p-value < 2.2e-16

ad.test(dt$tni_conc[dt$gender == "M"])


    Anderson-Darling normality test

data:  dt$tni_conc[dt$gender == "M"]
A = 41.122, p-value < 2.2e-16

# Kolmogorov-Smirnov test
lillie.test(dt$tni_conc[dt$gender == "F"])


    Lilliefors (Kolmogorov-Smirnov) normality test

data:  dt$tni_conc[dt$gender == "F"]
D = 0.30868, p-value < 2.2e-16

lillie.test(dt$tni_conc[dt$gender == "M"])


    Lilliefors (Kolmogorov-Smirnov) normality test

data:  dt$tni_conc[dt$gender == "M"]
D = 0.3578, p-value < 2.2e-16

# Pearson test
pearson.test(dt$tni_conc[dt$gender == "F"])


    Pearson chi-square normality test

data:  dt$tni_conc[dt$gender == "F"]
P = 540.41, p-value < 2.2e-16

pearson.test(dt$tni_conc[dt$gender == "M"])


    Pearson chi-square normality test

data:  dt$tni_conc[dt$gender == "M"]
P = 822.02, p-value < 2.2e-16

# Cramer-Vonmises test
cvm.test(dt$tni_conc[dt$gender == "F"])

Warning in cvm.test(dt$tni_conc[dt$gender == "F"]): p-value is smaller than
7.37e-10, cannot be computed more accurately


    Cramer-von Mises normality test

data:  dt$tni_conc[dt$gender == "F"]
W = 6.0311, p-value = 7.37e-10

cvm.test(dt$tni_conc[dt$gender == "M"])

Warning in cvm.test(dt$tni_conc[dt$gender == "M"]): p-value is smaller than
7.37e-10, cannot be computed more accurately


    Cramer-von Mises normality test

data:  dt$tni_conc[dt$gender == "M"]
W = 8.2575, p-value = 7.37e-10

Variance test

levene.test(dt$tni_conc, group = dt$gender)


    Modified robust Brown-Forsythe Levene-type test based on the absolute
    deviations from the median

data:  dt$tni_conc
Test Statistic = 3.5679, p-value = 0.0597

Wilcox or Mann-Whitney test

# Wilcox test
wilcox.test(dt$tni_conc ~ dt$gender,
            alternative = "two.sided",
            paired = FALSE, conf.int = TRUE,
            conf.level = 0.95)


    Wilcoxon rank sum test with continuity correction

data:  dt$tni_conc by dt$gender
W = 9690, p-value = 1.324e-12
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -0.4386174 -0.2453333
sample estimates:
difference in location 
             -0.336053

# Storing the results from a Wilcox test
test <- wilcox.test(dt$tni_conc ~ dt$gender,
            alternative = "two.sided",
            paired = FALSE, conf.int = TRUE,
            conf.level = 0.95)
test


    Wilcoxon rank sum test with continuity correction

data:  dt$tni_conc by dt$gender
W = 9690, p-value = 1.324e-12
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -0.4386174 -0.2453333
sample estimates:
difference in location 
             -0.336053

Brunner-Munzel test

# Brunner-Munzel test



test <- npar.t.test(tni_conc ~ gender,
                    data = dt,
                    method = "t.app",
                    alternative = "two.sided",
                    info = FALSE)
summary(test)


 #-----Nonparametric Test Procedures and Confidence Intervals for relative  effects-----# 
 
 - Alternative Hypothesis:  True relative effect p is less or equal than 1/2 
 - Confidence level: 95 % 
 - Method = Brunner - Munzel - T - Approx with 318.265 DF 
 #---------------------------Interpretation---------------------------------------------# 
 p(a,b) > 1/2 : b tends to be larger than a 
 #--------------------------------------------------------------------------------------# 
 
 #----Data Info-------------------------------------------------------------------------# 
  Sample Size
F      F  182
M      M  186

 #----Analysis--------------------------------------------------------------------------# 
  Effect Estimator Lower Upper     T p.Value
1 p(F,M)     0.714  0.66 0.767 7.835       0

Adding p-values to a box-plot

#install.packages("ggsignif")
library(ggsignif)

test <- t.test(tni_conc ~ gender, 
               alternative = "two.sided", var.equal = FALSE,
               paired = FALSE, data = dt)


# Adding p-values to plots
ls(test)

 [1] "alternative" "conf.int"    "data.name"   "estimate"    "method"     
 [6] "null.value"  "p.value"     "parameter"   "statistic"   "stderr"

test$p.value

[1] 0.01085427

ggplot(dt, aes(y = tni_conc, x = gender, fill = gender)) +
  stat_boxplot(geom = "errorbar", width = 0.1) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = 0.4)) +
  theme_bw() +
  scale_fill_brewer("Gender", palette = "Set1") +
  scale_y_continuous(trans = 'log10', limits = c(0.05, 120)) +
  labs(x = "Gender", y = "Troponin I (log 10 scale)") +
  theme(legend.position = "none")  +
  geom_signif(y_position = 2,
              xmin = 1,
              xmax = 2,
              annotation =
                ifelse(test$p.value < 0.001,
                       "< 0.001",
                       round(test$p.value, 3)),
              tip_length = 0.05,
              color = "black")

Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Removed 1 rows containing non-finite values (stat_boxplot).

Warning: Removed 1 rows containing non-finite values (stat_signif).

Warning: Removed 1 rows containing missing values (geom_point).