IntroSlides.qmd

---
src: plugin/math/math.js
title: "IntroSlides with focus on statistics"
format: 
  revealjs:
    # theme: league
    slide-number: c/t
    width: 1200
    height: 900
    transition: slide
    cap-location: bottom
    fig-align: center
    fig-format: png 
    fig-dpi: 300
    center: false
    chalkboard: true
---

## Basic steps in workflow

1.  *Define environment*
2.  *Import*
3.  *Transform*
4.  *Explore (general/outlier/distribution) (go back to 3?)*
5.  *Classify scale level / distribution (based on 3/4)*
6.  *Describe*
7.  *Test / Model (may include step 6)*
8.  *Report*

## Define environment

-   *Activate packages to use: library() / pacman::p_load()*
-   *ggplot theme: theme_set() / theme_update()*
-   *flextable settings: set_flextable_defaults()*
-   *knitr::opts_chunk\$set()*

```{r, echo=TRUE}
if(!requireNamespace("pacman")){install.packages("pacman")}
pacman::p_load(conflicted,tidyverse,wrappedtools, readxl,car, flextable,
               ggbeeswarm, ggsignif, ggridges, patchwork,ggrepel, easystats)

conflicts_prefer(dplyr::filter,dplyr::select)
theme_set(theme_light(base_size = 20))
gdtools::register_gfont('Roboto')# Mono')
set_flextable_defaults(
  theme_fun = theme_zebra, font.size = 18, font.family = 'Roboto',
  table.layout = 'autofit',
  padding.bottom = .2, padding.top = .2, padding.left = 2, padding.right = 2)

knitr::opts_chunk$set(message = FALSE, warning = FALSE, comment = NA, echo = TRUE)
```

## Import

-   *read_xlsx() / read_csv() / read_csv2() / read_spss() / read_sas() / DBI::xx() / dbplyr::xx()*
-   options related to separators, number formats, ranges etc.
-   *rename() / rename_with()*

```{r}
rawdata <- read_excel('Data/DOC-20230130-WA0000_.xlsx',
                      sheet = 1,col_names = TRUE)
```

## Glimpse at data: Find the problems?

```{r}
head(rawdata,n = 15) |> flextable()|> 
  theme_zebra(even_body = 'aquamarine',odd_body = 'antiquewhite')
```

## Rename

```{r}
colnames(rawdata)
rawdata <- rawdata |> 
  rename(Region=`CODE OF SAMPLE`) |> 
  rename_with(.fn = ~str_replace_all(
    .,
    c("GTH"="GHT", 'AL.+UM'= 'Cup', 
      'C[UA]P' = 'Cup','\\(\\w+.*\\)'='',
      'Cup Cup'='Cup','  '=' ')) |> 
      str_to_title() |> str_trim())
cn() 
```

## Transform

-   *Change or create columns with mutate() / mutate(across())*
-   *e.g. for log-transformation, creation of factors, text recoding*

```{r}
rawdata <- rawdata |> 
  mutate(
    `Code Of Cup`=case_when(
      is.na(`Code Of Cup`)~paste("sample",row_number()),
      .default = `Code Of Cup`),
    `Weight Of Sample After Drying`=`Weight Of Cup + Sample After Drying`-
      `Weight Of Empty Cup`,
    `Dry Content (%)`=`Weight Of Sample After Drying`*100/
      `Weight Of Sample Before Drying`,
    `Moisture Content (%)`=100-`Dry Content (%)`)
```

```{r, echo=FALSE}
head(rawdata,5)|> flextable()|> 
  theme_zebra(even_body = 'aquamarine',odd_body = 'antiquewhite')
```

## Explore / group variables

### *Explore (general/outlier/distribution)*

-   *ggplot()+geom_boxplot() / geom_beeswarm() / geom_density()*
-   *ks.test() / ksnormal() / shapiro.test()*

### *Classify scale level / distribution*

-   *gaussvars / ordvars / factvars, possibly more...*
-   *Store variables accordingly, e.g. ColSeeker()*

## *Explore: Outlier*

```{r, echo=TRUE}
p1 <- ggplot(data = rawdata, aes(x = `Region`, y = `Weight Of Empty Cup`))+
  geom_boxplot(outlier.alpha = 0) + #hide outliers, beeswarm will plot them
  geom_beeswarm(alpha=.5, color="darkorange")
p2 <- ggplot(data = rawdata, aes(x = `Region`, y = `Weight Of Sample Before Drying`))+
  geom_boxplot(outlier.alpha = 0) + 
  geom_beeswarm(alpha=.5, color="darkorange") 
p2 <- label_outliers(p2, labelvar = "Code Of Cup")
p1|p2

```

## Handle outliers?

Removal is the worst strategy possible, correct errors, think about distributions, winsorize...

![](Figures/outlier.jpg)

## *Explore: Unexpecteds*

```{r}
ggplot(data = rawdata,
       aes(x = `Region`,
           y = `Moisture Content (%)`))+
  geom_boxplot(outlier.alpha = 0) + 
  geom_beeswarm(alpha=.5, size = 2.5, color="darkorange")
```

## Transform Subregions?

```{r}
rawdata <- mutate(rawdata, Subregion = case_when(
  `Region`=='D' ~ 'D',
  `Region`=='A' & `Moisture Content (%)` > 90 ~ 'A high',
  `Region`=='A' & `Moisture Content (%)` <= 90 ~ 'A low') |>
    factor())
# Test for Region A is redundant here, but more verbose.

ggplot(data = rawdata,aes(x = `Moisture Content (%)`,y=Subregion))+
  geom_density_ridges(alpha=.5, scale=1)
```

## *Explore: Normal distribution 1*

-   Gaussian Normal distribution is required for many statistical procedures

-   Common tests are graphical exploration, Shapiro-Wilk-test and Kolmogorov-Smirnov-test

```{r, echo=TRUE}
p_normal <- 
  shapiro.test(x = rawdata$`Moisture Content (%)`)
p_normal

```

## 

```{r}

ggplot(data = rawdata,aes(x = `Moisture Content (%)`,fill=`Subregion`))+
  geom_density(linewidth=3, alpha=.5)+
  labs(title = paste('p (Shapiro) global',
                formatP(pIn = p_normal$p.value, pretext = T)),
       subtitle = rawdata |> group_by(Subregion) |> 
  summarize(pNormal=shapiro.test(`Moisture Content (%)`)$p.value |> formatP()) |> 
  unite('p(Normal)',sep = ': p=') |> pull(1) |> paste(collapse = '; '))
```

## *Explore: Normal distribution 2*

```{r, echo=TRUE}
rawdata |> 
  group_by(`Subregion`) |> 
  summarize(across(.cols = where(is.numeric), 
                   .fns = ~ksnormal(.x) |> #computing p-value
                     formatP(mark = T))) |> #formatting p-value
  pivot_longer(cols = -1, names_to = 'Variable', 
               values_to = 'pKS') |> #intermediate, all p-values in 1 column
  pivot_wider(names_from = `Subregion`, 
              values_from = pKS) #spreading across subregions
```

## *Explore: Group variables by type/distribution*

Scale level determines what statistics are appropriate

Typical scale levels are

-   nominal/categorical/factorial/qualitative: just different groups\
    (species, eye color, genotype, treatment)

-   ordered categories: few groups with inherent order\
    (quality bad\<medium\<good, pain between 0 and 10)

-   ordinal measures: many different values, natural order, no distribution assumption (satisfaction on a scale from 0 to 100)

-   measures following a Normal distribution

-   possibly measures from other known distributions\
    (beta, log-normal, poisson...), often treated as ordinal

## *Make type decision obvious/reproducible*

```{r, echo=TRUE}
gaussvars <- ColSeeker(data=rawdata,namepattern = c('Weight','Content'))
gaussvars
ordvars <- ColSeeker(namepattern='Weight.+Sample', exclude = 'After')
ordvars$names
factvars <- ColSeeker(namepattern='region',casesensitive = FALSE)
factvars$bticked
```

## Model

### *Describe*

-   *mean() / sd() / meansd()*
-   *median() / quantile() / median_quart()*
-   *table() / prop.table() / cat_desc_stats()*

### *Test*

-   *t.test() / lm()+\[Aa\]nova() / compare2numvars()*
-   *wilcox.test()*
-   *fisher.test() / glm(family=binomial)*

## *Model: Describe*

|  |  |
|------------------------------------|------------------------------------|
| **Sample size n:** per variable, if there are NAs |  |
| **Mean**: central tendency, the expected *typical* value | $$\frac{\sum{x}}{n}$$ |
| **Variance**: measure for variability/heterogeneity of data | $$\frac{\sum{(x-mean)^2}}{n-1}$$ |
| **Standard deviation SD**: the *typical* weighted deviation from the mean | $$\sqrt{Var}$$ |

## 

|  |  |
|------------------------------------|------------------------------------|
| **Standard error of the mean SEM**: how reliable is the mean *estimate*, what would be the expected SD of means from repeated experiments? | $$\frac{SD}{\sqrt{n}}$$ |
| **Median**: Split between lower/upper 50% of data |  |
| **Quartiles**: Split at 25%/50%/75% of data (more general: **Quantiles**, e.g.**Percentiles**), used in boxplot | various computational approaches |

## 

```{r, echo=TRUE}
desc_gauss <- rawdata |> 
  summarize(across(.cols = gaussvars$names,
                   .fns = meansd))
desc_gauss
desc_ord <- rawdata |> 
  summarize(across(ordvars$names,.fns=~median_quart(.x,roundDig = 3))) |> 
  pivot_longer(everything(),
               names_to = 'Measure', values_to = 'Median[1Q/3Q]')
desc_ord
```

## Descriptive Stats should match distribution and data

![](Figures/MeanMedian.jpg){fig-align="center"}

## *Model: Test*

Tests require hypotheses

![](Figures/guestimate.gif){fig-align="center"}

## Null hypothesis ?

-   Working hypothesis: This is what you expect!\
    E.g. treatment is lowering blood pressure more than placebo, transgenic animals become obese, bio reactor A is more efficient than B, concentration of substance is correlated with speed of reaction ...

-   Null hypothesis: This is what you test!\
    No difference / relation, BP under therapy = BP under placebo

## 4 possibilities:

-   Null hypothesis correct, test false positive (case A): alpha-error

-   Null hypothesis correct, test correct negative (case B)

-   Null hypothesis false, test false negative (case C): beta-error

-   Null hypothesis false, test correct positive (case D)

**Significance**: NOT probability of case A, but probability of your data given the NULL hypothesis, calculated from your data, conventionally \<0.05

**Power**: Probability of case D, *estimated* based on assumptions about effects and sample size, *calculation* would require knowledge of true difference, conventionally set at 0.80

## Test functions

### t-test / Wilcoxon-test (aka Mann-Whitney U-test)

![](Figures/MeanDiff.png){fig-align="center" width="800"}

## t-test

-   Assumptions: Continuous data with Normal distribution

-   1 or 2 (independent or dependent) samples with/without equal variances

-   how big is the mean difference relative to uncertainty?\
    t = (mean~1~ - mean~2~)/SEM

-   t follows a t-distribution, allows estimation of probability of t under the NULL hypothesis

### Wilcoxon-test

-   nonparametric, no distribution is assumed

-   based on rank-transformed data

-   insensitive to extreme values

## Test examples: *single variables*

```{r, echo=TRUE}
#t-Test with test for equal variances
t_out <- t.test(formula=`Moisture Content (%)`~`Region`, data=rawdata,
                var.equal=var.test(
                  formula=`Moisture Content (%)`~`Region`, 
                  data=rawdata)$p.value>.05)
t_out
#Wilcoxon-Test
wilcox.test(`Moisture Content (%)`~`Region`, 
            data = rawdata)
```

## From t to p

```{r, echo=FALSE}
tdata <- tibble(t=seq(-3,3,.01),y=dt(t,df = t_out$parameter))
plottmp <- ggplot(tdata,aes(t,y))+geom_line()
d <- ggplot_build(plottmp)$data[[1]]
plottmp + geom_area(data = d %>% filter(x>=abs(t_out$statistic)),aes(x,y),
                    fill='orangered',alpha=.5)+
  geom_area(data = d %>% filter(x<=-abs(t_out$statistic)),aes(x,y),
            fill='orangered',alpha=.5)+
  geom_vline(xintercept = t_out$statistic, linetype=3)+
  scale_x_continuous(breaks=seq(-10,10,.5))+
  ylab('probability density')+
  ggtitle(paste0('from t-test: t = ',roundR(t_out$statistic), 
                 ', p ',formatP(t_out$p.value, pretext = T)))+
  annotate('label', x=t_out$statistic,y=.0,
           label=paste('area = ', 
                       roundR(pt(abs(t_out$statistic),t_out$parameter,lower.tail = F))), 
           hjust=-0.1, vjust=-0.1)

```

## *Model: Test 2 / multiple variables*

```{r, echo=TRUE}
test_gauss <- compare2numvars(data = rawdata,
                              dep_vars = gaussvars$names,
                              indep_var = 'Region',
                              gaussian = TRUE,
                              round_p = 5)
test_gauss |> flextable()|> 
  theme_zebra(even_body = 'aquamarine',odd_body = 'antiquewhite')

```

## 

```{r}
test_ord <- compare2numvars(data = rawdata,
                            dep_vars = ordvars$names,
                            indep_var = 'Region',
                            gaussian = FALSE,round_desc = 3)
test_ord |> flextable() |> 
  theme_zebra(even_body = 'aquamarine',odd_body = 'antiquewhite')
```

## Show results

```{r, echo=FALSE, include=FALSE}
t_out <- t.test(formula=`Weight Of Sample Before Drying`~`Region`, data=rawdata,
                var.equal=var.test(
                  formula=`Weight Of Sample Before Drying`~`Region`, 
                  data=rawdata)$p.value>.05)
```

```{r}
ggplot(rawdata, aes(x = `Region`,y = `Weight Of Sample Before Drying`))+
  geom_boxplot(outlier.alpha = 0)+
  geom_beeswarm(alpha=.7, size=2.5,cex = 1.5, color="darkorange")+
  annotate(geom = 'label',x=2,y=4.2,
           label='extreme value -> high SD -> high SEM ->low t-statistic->high p-value',
           hjust=0.8,vjust=0)+
  geom_signif(comparisons = list(c(1,2)),
              annotations = paste('p',formatP(t_out$p.value)))
```

## Re-thinking test decision?

```{r}
ggplot(rawdata, aes(x = `Region`,y = `Weight Of Sample Before Drying`))+
  geom_boxplot(outlier.alpha = 0)+
  geom_beeswarm(alpha=.7, size = 2.5,cex = 1.5, color="darkorange")+
  geom_signif(comparisons = list(c(1,2)),test = wilcox.test)
```

## *Model: linear models 1 / univariable*

```{r}
plottmp <- ggplot(rawdata,aes(Subregion,`Moisture Content (%)`))+
  geom_violin(draw_quantiles = c(.25,.5,.75))+
  geom_beeswarm(color="darkorange")+
  
  ggtitle('Are all Subregions equal?')
print(plottmp)
```

## ANOVA: build model

```{r}
rawdata |> group_by(Subregion) |> 
  summarize(MeanMoisture=mean(`Moisture Content (%)`) |> roundR(4)) |> 
  pivot_wider(names_from = Subregion,values_from = MeanMoisture) |> 
  rename_with(~paste('Mean moisture %\n',.x)) |> flextable()|>
  theme_zebra(even_body = 'aquamarine',odd_body = 'antiquewhite')
lm1<- lm(`Moisture Content (%)`~Subregion, data=rawdata)
lm1
```

## ANOVA: get p-values

```{r}
anova(lm1) |> broom::tidy() |> flextable()|> 
  theme_zebra(even_body = 'aquamarine',odd_body = 'antiquewhite')
#post-hoc
(posthoc_out <- pairwise.t.test(x = rawdata$`Moisture Content (%)`, 
                                g = rawdata$Subregion, 
                                p.adjust.method = 'fdr')$p.value |> 
    formatP(ndigits = 5))
```

<!-- ## *Model: linear models 2 / multivariable* -->

```{r, echo=FALSE, eval=FALSE, include=FALSE}
lm2<- lm(sysRR~(Sex+Agegroup)*Treatment, 
         data=faketrial)
lm2
Anova(lm2,type = 3)
```

## *Visualize ANOVA*

```{r}
ggplot(rawdata,aes(Subregion,`Moisture Content (%)`))+
  geom_violin(draw_quantiles = c(.25,.5,.75))+
  geom_beeswarm(color="darkorange")+
  geom_signif(comparisons = list(c(1,2),c(1,3),c(2,3)),
              annotations = c(posthoc_out[,1], posthoc_out[2,2]),
              step_increase = .25)+
  scale_y_continuous(expand = expansion(mult = .1))+
  ggtitle('All Subregions are unequal')
```

## *Analyze more than 1 outcome*

```{r}
test_out <- compare_n_numvars(.data=rawdata,
                              dep_vars=gaussvars$names,
                              indep_var='Subregion',
                              gaussian=TRUE)
test_out$results |> 
  select(Variable,contains("fn"),multivar_p) |> 
  rename_with(~str_remove(.x," fn")) |> 
  flextable()|> 
  theme_zebra(even_body = 'aquamarine',odd_body = 'antiquewhite') |> 
  add_footer_lines(
    values='b indicates difference from group 2, c indicates difference from group 3')

```

## Regression: Scatterplot

```{r, fig.width=15, fig.height=10, echo=FALSE}
p1 <- ggplot(rawdata,aes(`Weight Of Empty Cup`,`Weight Of Sample Before Drying`))+
  geom_point(color = "darkorange", size = 2, alpha=.4) + 
  geom_smooth()
p2 <- rawdata |> filter(`Weight Of Sample Before Drying`>4.5) |> 
  ggplot(aes(`Weight Of Empty Cup`,`Weight Of Sample Before Drying`))+
  geom_point(color = "darkorange", size = 2, alpha=.4) + 
  geom_smooth()+
  geom_smooth(method='lm')
p3 <- rawdata |> filter(`Weight Of Sample Before Drying`>4.5) |> 
  ggplot(aes(`Weight Of Empty Cup`,`Moisture Content (%)`))+
  geom_point(color = "darkorange", size = 2, alpha=.4) + 
  geom_smooth(method='lm')
p1+p2+p3
```

## Regression: Underlying *mechanics*

![](Figures/RegressionSprings.png){fig-align="center"}

## Regression: Statistics

```{r}
lm_out0 <- lm(`Weight Of Sample Before Drying` ~ `Weight Of Empty Cup`, 
             data=rawdata)
lm_out0

# filtering outlier
lm_out <- lm(`Weight Of Sample Before Drying` ~ `Weight Of Empty Cup`, 
             data=rawdata |> filter(`Weight Of Sample Before Drying`>4.5))
lm_out

```

## *Regression: Visualize*

```{r, echo=FALSE, fig.width=15}
ggplot(rawdata,aes(`Weight Of Empty Cup`,
                   `Weight Of Sample Before Drying`))+
  geom_point(aes(color=`Weight Of Sample Before Drying`>4.5),alpha=.5)+
  scale_color_manual(values=c('darkred','darkorange'))+
  guides(color='none')+
  geom_smooth(data = rawdata |> filter(`Weight Of Sample Before Drying`>4.5),
              method='lm', se=FALSE)+
  geom_abline(slope=lm_out$coefficients[2],
              intercept = lm_out$coefficients[1],
              linetype=3, color='darkred')+
  scale_x_continuous(limits=c(-1,NA))+
  scale_y_continuous(limits=c(NA,lm_out$coefficients[1]))+
  annotate('label',x=0,y=lm_out$coefficients[1], size=7,
           label=paste('Intercept:',roundR(lm_out$coefficients[1],4)),hjust=-.1)+
  geom_segment(x=0,xend=1,
               y=lm_out$coefficients[1]+lm_out$coefficients[2],
               yend=lm_out$coefficients[1]+lm_out$coefficients[2],
               linetype=2)+
  geom_segment(x=0,xend=0,y=lm_out$coefficients[1],
               yend=lm_out$coefficients[1]+lm_out$coefficients[2],
               linetype=2,
               arrow = arrow(length = unit(0.2, "inches"), 
                             type='closed', ends='both'))+
  annotate('label',x=0,y=lm_out$coefficients[1]+lm_out$coefficients[2]/2, size=7,
           label=paste('Slope:',roundR(lm_out$coefficients[2],4)),hjust=1.1)
```

## Regression: Significance

```{r}
anova(lm_out) |> broom::tidy()
model_parameters(lm_out)
```

## *Report*

-   *RMarkdown and quarto are powerful tools to create reports and presentations*

-   Export figures: ggsave() / png() / pdf()

-   Export tables: write_xlsx()

-   *Package flextable provides nice features for table formatting*

## Flextable example

```{r, echo=TRUE,size=3}
test_ord |> select(-desc_all) |> rename_with(~str_remove(.,'Code Of ')) |> 
  flextable() |> 
  theme_zebra(even_body = 'aquamarine',odd_body = 'antiquewhite')|> 
  italic(~p<=0.05,j = 1) |> bg(~p<=0.05,j = 4,bg = 'yellow') |> 
  set_caption('Treatment effects, measures following a normal distribution') |> 
  add_footer_lines('Significance level is set at 0.05') |> 
  fontsize(size = 12,part = 'footer') 

```

## Useful tools along the way

-   Pick columns / rows: select() / pull() / filter() / slice()

-   Change format of tibble wide \<--\> long (e.g. for repeated measures):\
    pivot_longer()/pivot_wider()

-   Regular expressions: str_replace() / str_detect() / str\_...

-   Merge text elements: paste() / str_glue()

-   Apply functions: purrr::map_xxx