Statistical Analysis of Biological Data

Link to test summary slide

• Summary slide for part 3 stat

Background

In this course we will explain the theory behind common statistical tests to compare two or more groups. We will illustrate this using biologically-motivated examples and R code that uses the “tidyverse” techniques of data manipulation and visualization.

To install the entire collection of tidyverse packages (which includes other useful data-related packages that we will not use today), you may wish to use the following:-

## May take some time!
install.packages("tidyverse")

Statistics has been a fundamental part of the R language from the beginning. A vast array of statistical tests and associated visualizations are supported. However, the functions to perform these predate the creation of the tidyverse and are used in a manner that might be somewhat counter-intuitive to someone that is familiar with the tidyverse. Therefore we will use a couple of add-on packages that have been created to perform statistics in a manner that is compatible with the tidyverse eco-system.

• rstatix: Statistical testing for the tidyverse
• ggpubr: ‘ggplot2’ Based Publication Ready Plots

The specific set of required packages is as follows:-

install.packages(c("dplyr",
                   "ggplot2",
                   "readr",
                   "readxl",
                   "rstatix",
                   "ggpubr",
                   "rmarkdown",
                   "tidyr",
                   "vcd")) 

Part I - Contingency tables

When working with categorical variables, we are usually interested in the frequencies of the different categories in our sample. To display data for two or more categorical variables, cross-tabulations, or contingency tables, are commonly used - with 2 x 2 tables being the simplest. We can then test whether there is an association between the row factor and the column factor by a chi-squared test or a Fisher’s exact test.

To demonstrate the analysis of contingency tables we will use a dataset provided with the vcd package. You will need to install this package using the install.packages R function.

install.packages("vcd")

## Load the packages we will need
library(dplyr)
library(ggplot2)
library(rstatix)
library(ggpubr)
library(vcd)

The data frame Arthritis should then be accessible which is described as:-

Data from Koch & Edwards (1988) from a double-blind clinical trial investigating a new treatment for rheumatoid arthritis.

#let's use the'Arthritis' dataset in the 'vcd' package 
Arthritis

The count function, included in dplyr can be used to give a tabulation of the values in any column in the data frame.

##one way table (count of one variable)
count(Arthritis,Sex)

We can quickly visualise these counts as a barplot with ggplot2. The geom_bar plot type will automatically count the number of observations that will be plotted on the y-axis.

ggplot(Arthritis, aes(x = Sex)) + geom_bar()

The function freq_table from rstatix is another way of making the counts and will also add the proportions as an extra column.

#to get proportion of males and females
freq_table(Arthritis,Sex) 

The count function can also compare two columns from a data frame if both columns are given as arguments.

count(Arthritis,Sex,Improved)

These can also be plotted, but this time we have to use geom_col to specify the values for the y-axis.

count(Arthritis,Sex,Improved) %>% 
ggplot(aes(x = Sex, fill = Improved, y = n)) + geom_col()

count(Arthritis,Sex,Improved) %>% 
ggplot(aes(x = Sex, fill = Improved, y = n)) + geom_col(position = "dodge")

The freq_table function can be used to calculate the proportions. However, we need to pay attention to the order in which the variables are specified as this will dictate how the proportions are calculated

Compare the output of:-

freq_table(Arthritis,Sex, Improved)

to:-

freq_table(Arthritis,Improved, Sex)

Having explored our data, we can now perform statistical testing. The function chisq_test can be used to assess whether differences in proportions are significant or not. We actually don’t need to calculate the proportions; R will do this for us.

However, we need to re-format the data slightly into a wide table rather than the default long nature of a data frame in the tidyverse. We the pivot_wider function we create a two-by-two table that is typically used for a contingency analysis.

Arthritis %>% 
  count(Improved, Treatment) %>% 
  tidyr::pivot_wider(values_from = n,names_from = Improved)

##would also work with names_from Treatment

And now remove the Treatment column as the chisq_test function is only expecting numeric data. The results are presented in a tidyverse tibble that and we can interpret the test statistics and p-value from this table

Arthritis %>% 
  count(Improved, Treatment) %>% 
  tidyr::pivot_wider(values_from = n,names_from = Improved) %>% 
select(-Treatment)  %>% 
  chisq_test()

Note that the only line of code to perform the test is chisq_test. The preceding lines are used to manipulate our data into the correct format. If our data are already in a numeric table we might be able to use chisq_test directly.

The chi-squared test works by comparing the frequencies in each cell of our table to what we would expect by chance (i.e. if there was no significant association). If we want to see what the expected frequencies would be we can run expected_freq after chisq_test.

Arthritis %>% 
  count(Improved, Treatment) %>% 
  tidyr::pivot_wider(values_from = n,names_from = Improved) %>%
select(-Treatment)  %>% 
  chisq_test() %>% expected_freq()

     None     Some   Marked
[1,] 21.5 7.166667 14.33333
[2,] 20.5 6.833333 13.66667

However, the chisq_test function is not appropriate in all circumstances.

Arthritis %>% 
  count(Improved, Sex) %>% 
  tidyr::pivot_wider(values_from = n,names_from = Improved) %>% 
  select(-Sex)  %>% 
  chisq_test 

Warning: Chi-squared approximation may be incorrect

The Fisher test is recommended for tables with low numbers of observations (e.g. when more than 20% of cells have expected frequencies < 5)

Arthritis %>% 
  count(Improved, Sex) %>% 
  tidyr::pivot_wider(values_from = n,names_from = Improved) %>% 
 select(-Sex)  %>% 
  fisher_test

In this case we had some low expected frequencies so the Fisher test was more appropriate

Arthritis %>% 
  count(Improved, Sex) %>% 
  tidyr::pivot_wider(values_from = n,names_from = Improved) %>% 
  select(-Sex)  %>% 
  chisq_test %>% 
  expected_freq()

Warning: Chi-squared approximation may be incorrect

     None     Some    Marked
[1,] 29.5 9.833333 19.666667
[2,] 12.5 4.166667  8.333333

Exercise

1- Read the excel file called Ex Biostat P1.xlsx into R (see below for the required code). We recommend using make_clean_names function from rstatix to make sure the column names are “clean” (i.e. without spaces or other characters that could cause issues for R).

2- Use the counts function to make a cross-tabulation of Tumor grade against Gender

3- Determine the proportion of Grade III tumors within females

4- Use the appropriate test to check if the tumor grade depends on the gender

## the readxl package is required to read xls and xlsx files into R
## However, csv and tsv files are recommended to store data
tab <- readxl::read_xlsx("data/EX Biostat P1.xlsx") %>% 
  make_clean_names()

Part II - How to assess normality

We will read some example data to illustrate how one would test for a normally-distributed variable. This property is important as it influences which test we should use.

One of the best ways of displaying data is by using a graph. Graphs can make both simple and complex data easier to understand by making it easier to spot trends and patterns. We can use plots to view the distribution of our data (minimum, maximum, mid-point, spread etc) and to ensure that the values in our dataset seem realistic (e.g. no outliers). Many statistical tests rely on the assumption that the data are normally distributed.

The data for this section are to be found in the file normal_example.csv in the data folder. You will need to specify the file path accordingly.

library(readr)
df1 <- read_csv("data/normal_example.csv")

We can inspect the data in RStudio and discover that it consists of a tidy dataset with numeric values in a column called Values and a column Var to indicate a variable name (x)

View(df1)

Various graphical methods are available to assess the distribution. The first of which is a histogram. In this graph the data are split into “bins” and the value on the y axis corresponds to the number of observations in that bin. The user only has to specify the variable to be plotted, and function takes care of the binning. From this plot we can judge what the average value of the data is, and the spread.

Histograms can be made in ggplot2 by using the geom_hist function (or indeed the base hist function). Here we are going to make use of the gghistogram function from ggpubr as it is a bit more convenient.

gghistogram(df1,x="Value")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

When assessing the distribution of a variable, you might be tempted to plot the histogram and density on the same plot. gghistogram has the argument add_density, which should do the job.

gghistogram(df1,x="Value",add_density = TRUE)

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

However, this doesn’t work. If you check the y-axis limits on the histogram and density plots, you’ll notice they are on a different scale. Conveniently there is an argument in gghistogram that solves this issue.

gghistogram(df1,x="Value",
            add_density = TRUE,
            y="..density..")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

We can add the standard normal curve with the following code. To help distingush the standard normal we can change the colour this is plotted in.

gghistogram(df1,x="Value",add_density = TRUE,y="..density..") + stat_overlay_normal_density(col="steelblue",lwd=2,lty=2)

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

A box plot is an excellent way of displaying continuous data when you are interested in the spread of your data. The “box” of the box plot corresponds to the lower and upper quartiles of the respective observations and the bar within the box, the median. The whiskers of the box plot correspond to the distance between the lower/upper quartile and the smaller of: the smaller/largest measurement OR 1.5 times the inter quartile range. A disadvantage of the box plot is that you don’t see the exact data points. However, box plots are very useful in large datasets where plotting all of the data may give an unclear picture of the shape of your data.

Again, we use ggboxplot from ggpubr for convenience.

ggboxplot(df1, y = "Value")

A violin plot is sometimes instead of the boxplot to provide more information about the density.

ggviolin(df1, y = "Value")

Individual points can also be added with the jitter argument; avoiding over-plotting by adding random noise along the x-axis.

ggviolin(df1, y = "Value",add = "jitter")

Finally, we have a “qq-plot” which allows to compare the quantiles of our dataset against a theoretical normal distribution. If the majority of points lie on a diagonal line then the data are approximately normal. The ggqqplot function in ggpubr is the most convenient way of creating this plot.

ggqqplot(df1, x="Value")

These graphical methods are by far the easiest way to assess if a given dataset is normally-distributed. However, you do not necessarily need to generate and report all the plots for your data; just as many as you to inform your decision.

For “real-life” data, the results are unlikely to give a perfect plot, so some degree of judgement and prior experience with the data type are required. Indeed, it should be noted that the dataset visualised in the above plots was sampled from a normal distribution. Even then, the plots were not 100% convincing!

Tests for normality

Although their usage is contentious amongst statisticians, there are a few methods for testing whether variables are normally-distributed or not. If the p-value is sufficiently small from these methods then we conclude that the data are not normally distributed. However, some statisticians prefer to use graphical methods and their intuition about the data or prior knowledge of the data type (e.g. some measures are generally believed to be normally-distributed)

#shapiro test from rstatix
#p<0.05 ...difference between data and normality..data not normal
#p>0.05 ...no diff between data and normality ..data normally distributed
shapiro_test(df1,Value) 

Descriptive Statistics

When performing a statistical analysis it is common practice to report on the average and variability. Our decision about whether the data are normally-distributed will influence what measures we report. In rstatix there is a function called get_summary_stats that can calculate such summaries.

For a dataset that is normally-distributed, appropriate measures of the average and variability are the mean and standard deviation.

df1 %>% 
  get_summary_stats(type="common") 

df1 %>% 
  get_summary_stats(type = "common") %>% 
  select(mean, sd)

Exercise

1- Read the excel file called Ex Biostat P2.xlsx into R

2- Decide if the age and hospitalization days are normally distributed. Vote for your findings on wooclap

3- Calculate the appropriate descriptive statistics [mean and SD, or median and IQR] for each variable

Part III - Significance tests for continuous variables

In this part we will show how to perform tests to compare 1, 2 (or more) continuous variables. The dataset, provided by MASH at The University of Sheffield, describes individuals that have been following different diets and their age and gender. The main goal of interest is to determine which of three competing diet regimes results in the greatest weight loss. However, we can use the dataset to demonstrate other types of test.

diet <- read_csv("data/diet.csv")
diet

One-sample test

The first hypothesis we will test is whether the people in the study are overweight or not. This first involves some manipulation of the table to calculate an extra variable; the Body Mass Index (BMI). We will test if people in our study are overweight, where overweight is defined as having a BMI over 25.

\(BMI = weight / height^2\)

(where the weight is measured in kg, and the height in metres)

Exercise (short) - Add a new variable to the data frame for the BMI of each person + you might want to do this in multiple steps using the %>% notation

diet <- diet %>% 
  mutate(BMI =initial.weight / (height/100)^2)

We can now test our new variable for normality using the plots and tests from earlier, although we will not show all the plots here.

gghistogram(diet, x = "BMI", ,add_density = TRUE,y="..density..")+ stat_overlay_normal_density(col="steelblue")

The one-sample t-test is implemented in the function t_test (as are the various types of t-test that we will see). The variables to be used in the test are defined using R’s formula ~ syntax.

Y ~ X

Where Y is a numeric variable and X is a categorical variable indicating the groups that each particular value of Y belongs to.

In the case of a one-sample test we are testing one numeric variable against a known or population mean (mu). As we don’t have a groups to compare, this is written as:-

Y ~ 1

diet %>% 
  t_test(BMI ~1)

The various statistical tests in rstatix all have the _test suffix (e.g. t_test, chisq_test, shapiro_test). The naming of these functions is presumably chosen to be consisent with the base R implementations which are called t.test chisq.test etc. However, the .test functions are older and do not work well with tidy data. Make sure you are using the correct function.

We get a hugely significant result! However, if we look at the description for t_test it is testing against a population mean of \(0\). It is no surprise that we get a significant result! By changing the mu argument we can perform a test to see if the people in the study are overweight to begin with (using 25 as the population or known mean).

diet %>% 
  t_test(BMI ~1,mu = 25,alternative = "greater")

## 25 is the cutoff for overweight

Two-sample tests

A two-sample t-test should be used if you want to compare the measurements of two populations. There are two types of two-sample t-test: independent (unpaired) and paired (dependent). To make the correct choice, you need to understand your underlying data.

• An independent two-sample t-test is used when the two samples are independent of each other, e.g. comparing the mean response of two groups of patients on treatment vs. control in a clinical trial.

• As the name suggests,a paired two-sample t-test is used when the two samples are paired, e.g. comparing the mean blood pressure of patients before and after treatment (two measurements per patient).

Back to our dataset, we might wonder if there is actually any effect due to diet so we will compare the intial and final weights.

Re-formating the dataset for tidy analysis

The dataset in it’s current form is not suitable for analysis with tidy methods. Before proceeding we need to re-format the data into two columns.

• One column containing numeric variables (each value of weight)
• An indicator (or categorical) variable to denote the group each numeric observation belongs to (initial or final weight)

The pivot_longer function supports this transformation. As we want to use all columns in our existing dataset we have to use the everything() shortcut to select the columns.

diet_long <- diet %>% 
  select(contains("weight")) %>% 
  tidyr::pivot_longer(everything(),names_to = "time_point", values_to = "weight") 

diet_long

So now weight contains all our observations of weight and time_point indicates when the weights were measured. The gghistogram function we introduced earlier is able to visualise data in this form, and we can use the facet.by argument to produce separate plots for each type of weight measurement.

diet_long %>% 
  gghistogram(x = "weight", facet.by = "time_point",add_density = TRUE, y = "..density..") + stat_overlay_normal_density(col="red")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

We can run the shapiro test on the two variables:-

diet_long %>%
  group_by(time_point) %>% 
  shapiro_test(weight)

The t_test function requires us to create a formula as before. This time we are doing a two-sided test, so the formula is Y ~ X where ‘Y’ is the numeric variable and ‘X’ is the categorical/groups variable.

diet_long %>%
  t_test(weight ~ time_point)

The p-value is significant and shows that overall the weights of individuals is different before and after diet. However, this test is not specifically testing for a decrease after the diet (which we would really hope to be the case). By adding an extra argument alternative=less we get a different result

diet_long %>%
  t_test(weight ~ time_point, alternative = "less")

There is also extra information that we could employ; namely that the measurements of weight are made on the same person before and after dieting. This is a classic example of when to apply a paired t-test. Again, we do not need to use a different function to perform the test; only add an argument paired=TRUE to t_test.

diet_long %>%
  t_test(weight ~ time_point, alternative = "less",paired=TRUE)

An intuitive way of comparing the distributions and showing the test result on the plot is using a boxplot. For this we will revert to ggplot2, but we are to include a p-value using the stat_compare_means function (from ggpubr.

diet_long %>% 
  ggplot(aes(x = time_point, y = weight)) + geom_boxplot() + geom_jitter(width=0.1) + 
  stat_compare_means(method = "t.test",paired=TRUE, method.args = list(alternative = "less"))

?stat_compare_means

A convenient plot in ggpubr will allow us to visualise the paired differences

ggpaired(diet, cond1="initial.weight", cond2="final.weight",line.color = "grey") + stat_compare_means(paired=TRUE,method = "t.test")

The Independent t test, with two independent groups

Lets consider that we want to compare whether males or females lost more weight during the trial. So, let’s create a new variable called weight loss

diet <- diet %>% 
mutate(wt.loss = initial.weight - final.weight)

Here we have two groups (males and females), and these can be treated as independent variables as each group has different participants than the other group.

The null hypothesis for such a test would be that the weight loss is the same between groups male and female. We seek to evidence to reject this hypothesis by calculating a test statistic.

Remember, the t-test assumes normal distribution and equal variance

Firstly, we have to check for normality:-

diet %>% 
ggdensity("wt.loss",color = "gender") + stat_overlay_normal_density()

diet %>%  
  group_by(gender) %>% 
  shapiro_test(wt.loss)

Then we check for equal variances using the Levene test

diet %>%
  levene_test(wt.loss~ gender)

Warning: group coerced to factor.

#if p value is significant --> difference in variance (unequal variance) -->  use Welch t-test (R default) 
#if p value is not significant --> no difference in variance (equal variance) --> use t-test (add argument var.equal =TRUE)

We conclude that the variances are approximately the same and both variables are normally-distributed

We can use the t_test function to perform an independent test. As before, the formula argument to t_test is the R formula notation for the test being performed. However, this time we do not need the pivot_longer function as we already the gender and weight.loss variables in our data frame

The t_test function allows various type of test to be performed by changing the appropriate arguments (see the help for t_test for details (?t_test)). For instance, we can tell the test that we believe our variances are equal or not.

diet %>% 
  t_test(wt.loss ~ gender,var.equal = TRUE) 

NA

We can see that the t-statistic we observe is consistent with the null hypothesis, that the weight loss in males and females is the same. That is, the probability of observing a t-statistic of 0.2 or more, or -0.2 or less, is quite high.

This is not a significant result (p>0.05), so there is no evidence of a difference in weight loss between males and females

diet %>% 
  ggplot(aes(x = gender, y=wt.loss)) + geom_boxplot() + stat_compare_means(method="t.test")

Non- Parametric alternatives (e.g. the Wilcoxon test)

Being able to use the t_test relies on the your data being normally-distributed. If we do not sufficient confidence in this assumption, there are different statistical tests that can be applied. Rather than calculating and comparing the means and variances of different groups they are rank-based methods. However, they still come with a set of assumptions and involve the generation of test statistics and p-values.

Independent samples = Wilcoxon rank sum test (Mann Whitney U test)

This test has many different names including the Wilcoxon, Wilcoxon two sample test, Wilcoxon-Mann-Whitney, Wilcoxon rank sum and the Mann-Whitney-U test. However, this test should not be confused with the Wilcoxon signed rank test (which is used for paired tests). To avoid confusion this test is usually referred to as the Mann-Whitney U test, which is used when the dependent variable to be examined is continuous but the assumptions for parametric tests are violated.

Fortunately, the rstatix developers have made the function to do a Wilcox-test similar to doing a t_test. The difficulty is in choosing the correct test to apply - which R will not advise you on.

Let’s go back to our example of comparing weight loss between groups male and female. The equivalent non-parametric version of the test we performed before is:-

diet %>% 
  wilcox_test(wt.loss ~ gender)

The wilcox_test is flexible in much the same way that t_test is. We can switch to applying a paired test by adding the argument paired=TRUE.

diet_long %>% 
  wilcox_test(weight ~ time_point, paired=TRUE)

Compare between more than two groups

Parametric (ANOVA)

The two-sample t-test is useful when we have just two groups of continuous data to compare. When we want to compare more than two groups, a one-way ANOVA can be used to simultaneously compare all groups, rather than carrying out several individual two-sample t-tests. The main advantage of doing this is that it reduces the number of tests being carried out, meaning that the type I error rate (the probability of seeing a significant result just by chance) does not become inflated.

In order to justify if an ANOVA test is appropriate we have to test for normality.

#by histogram
diet %>%   
gghistogram(x = "wt.loss", facet.by = "diet.type", add_density = TRUE, y = "..density..") + stat_overlay_normal_density(col="red")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

#by Q-Q plot
diet %>% 
  ggqqplot(x = "wt.loss", facet.by = "diet.type")

#by shapiro test
group_by(diet, diet.type) %>% 
  shapiro_test(wt.loss)

A one-way ANOVA compares group means by partitioning the variation in the data into between group variance and within group variance. Like the other statistical tests we have encountered, the functions in R do the hard work of calculating the statistics. The anova_test function is a tidy version of the ANOVA test.

diet %>% 
  anova_test(wt.loss ~ diet.type)

ANOVA Table (type II tests)

     Effect DFn DFd     F     p p<.05   ges
1 diet.type   2  73 5.383 0.007     * 0.129

When the test provides a significant result (like above) it tells us that there is at least on difference in the groups. However, it does not tell us which group is different. For this, we can apply a “post-hoc test” such as the Tukey test. If anova_test did not produce a significant p-value, we wouldn’t proceed with this step

diet %>%  
  tukey_hsd(wt.loss ~ diet.type)

As we have seen previously, a standard method of presenting the differences between groups is to use the stat_compare_means function to automatically add p-values to a boxplot or violin plot

ggplot(diet, aes(x = diet.type, y = wt.loss)) + geom_violin() + geom_jitter(width=0.1) + stat_compare_means(method="anova")

However, in the case of more than two groups it will only show a single p-value from the ANOVA rather than individual comparisons. We can explicitly list particular contrasts we are interested in.

my_comparisons <- list( c("A", "B"), c("A", "C"), c("B", "C") )
ggplot(diet, aes(x = diet.type, y = wt.loss)) + geom_violin() + geom_jitter(width=0.1) + 
  stat_compare_means(method = "t.test",comparisons = my_comparisons)

Alternatively, we can manually-compute the p-values and add these to the plot.

stat_res <- diet %>% 
  tukey_hsd(wt.loss ~ diet.type)


ggplot(diet, aes(x = diet.type, y = wt.loss)) + geom_violin() + geom_jitter(width=0.1)+
  stat_pvalue_manual(stat_res, label = "p.adj",y.position = c(11, 13, 15))

Non-Parametric (Kruskal Wallis)

Data that do not meet the assumptions of ANOVA (e.g. normality) can be tested using a non-parametric alternative. The Kruskal-Wallis test is derived from the one-way ANOVA, but uses ranks rather than actual observations. It is also the extension of the Mann-Whitney U test to more than two groups.

diet %>% 
  kruskal_test(wt.loss ~ diet.type)

Like the one-way ANOVA this will only tell us that at least one group is different and not specifically which group(s). The post-hoc dunn.test is recommended which also performs multiple testing correction.

diet %>% 
  dunn_test(wt.loss ~ diet.type, p.adjust.method = "bonferroni")

At this point we could be about to recommend diet C to those that wish to lose weight. However, are there any other factors in the data that we should be considering? With ggplot2 we can quite easily visualise the effects of multiple factors on the data. Lets add both gender and diet type into the plot. It now appears that diet C is having an effect on males but not females.

ggplot(diet, aes(x = diet.type, y = wt.loss,fill=gender)) + geom_boxplot()

Two-way ANOVA

The formula notation allows us to specify an interaction between gender and diet type. In other words, we are looking to see if the effect of diet type is different for males and females. In R, the formula for an interaction is specified using a * between the variables that we are interested in assessing the interaction for.

diet %>% 
  anova_test(wt.loss ~ diet.type*gender)

ANOVA Table (type II tests)

            Effect DFn DFd     F     p p<.05      ges
1        diet.type   2  70 5.619 0.005     * 0.138000
2           gender   1  70 0.031 0.860       0.000448
3 diet.type:gender   2  70 3.153 0.049     * 0.083000

This tells us that an effect exists between diet type and gender, but like before we have to run a post-hoc test to discover more

diet %>% 
  tukey_hsd(wt.loss ~ diet.type*gender) %>% 
  filter(p.adj.signif != "ns")

Exercise

The excel file ‘RCC2’ contains data about the expression levels of some genes in patients with renal cell carcinoma. In your study, you put the following hypotheses. Please test those alternative hypotheses and state whether you will accept or reject each one.

1. Females have a higher level of E2F3 than males
2. ANXA expression levels vary between unilateral and bilateral tumors
3. Individuals with RCC grade II have different levels of E2F3 than those with grade III or IV
4. The mean value of miR499 decreases significantly after treatment
5. DFFA is higher in patients with grade IV tumors

Solutions

To be revealed during the workshop!

Exercise 1

## NEW

tab <- readxl::read_xlsx("data/EX Biostat P1.xlsx") %>% make_clean_names()
##2-Make a cross table showing the gender in the rows and tumor grade in the columns

count(tab, Gender, Tumor.grade)

##3-Define the percentage of Grade III tumors within females 
freq_table(tab, Gender, Tumor.grade)

##4-Use the appropriate test to check if the tumor grade depends on the gender 

count(tab, Gender, Tumor.grade) %>% 
    tidyr::pivot_wider(names_from = Gender,values_from = n) %>% 
  select(-Tumor.grade) %>% 
  chisq_test()

NA

Exercise 2

##1-Read the excel file called 'EX Biostat P2' into R 
Biostat2 <- readxl::read_xlsx("data/EX Biostat P2.xlsx")

gghistogram(Biostat2, x = "Age")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

gghistogram(Biostat2, x = "Age", add_density = TRUE, y = "..density..") + stat_overlay_normal_density(col="blue")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

ggqqplot(Biostat2, x = "Age")

Biostat2 %>% 
  shapiro_test(vars = "Age")

# Age is normally distributed

gghistogram(Biostat2, x = "hospitalization days")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

gghistogram(Biostat2, x = "hospitalization days", add_density = TRUE, y = "..density..") + stat_overlay_normal_density(col="blue")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

ggqqplot(Biostat2, x = "hospitalization days")

Biostat2 %>% 
  shapiro_test(vars = "hospitalization days")

# hospitalization days NOT normally distributed

get_summary_stats(Biostat2)

NA

Exercise 3

3.1

RCC2 <- readxl::read_xlsx("data/RCC2.xlsx")

gghistogram(RCC2, x = "E2F3", 
            y = "..density..",
            facet.by = "gender") + stat_overlay_normal_density(col="blue")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

RCC2 %>% 
  group_by(gender) %>% 
  shapiro_test(E2F3)

#both groups normally distributed
#parametric * 2 groups --> t-test

RCC2 %>% 
  group_by(gender) %>% 
  get_summary_stats("E2F3",type="common")

t_test(RCC2, E2F3 ~ gender,var.equal = TRUE)

#ANSWER=REJECT

ggplot(RCC2, aes(x = gender, y = E2F3)) + geom_violin() + geom_jitter(width = 0.1) + stat_compare_means(method = "t.test")

3.2

gghistogram(RCC2, x = "ANXA", facet.by = "Side")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

RCC2 %>% 
  group_by(Side) %>% 
  shapiro_test(ANXA)

#not normal
#non-parametric * 2 groups --> Wilcoxon (Mann Whitney U)
wilcox_test(RCC2, ANXA ~ Side)

#ANSWER=REJECT
ggplot(RCC2, aes(x = Side, y = ANXA)) + geom_violin() + geom_jitter() + stat_compare_means()

3.3

gghistogram(RCC2, x = "E2F3", facet.by = "Grade", y = "..density..") + stat_overlay_normal_density() 

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

#3 groups normal
#parametric * > 2 groups --> ANOVA

anova_test(RCC2, E2F3 ~ Grade)

ANOVA Table (type II tests)

  Effect DFn DFd     F     p p<.05  ges
1  Grade   2 111 5.492 0.005     * 0.09

#there is a difference but between which groups? --> post Hoc
tukey_hsd(RCC2, E2F3 ~ Grade)

#ANSWER= ACCEPT
my_comparisons <- list( c("II", "III"), c("II", "IV"), c("III", "IV") )

ggplot(RCC2, aes(x = Grade, y = E2F3)) + geom_violin() + geom_jitter(width=0.1) + stat_compare_means(method="t.test",comparisons = my_comparisons)

NA
NA

3.4

mir499_data <- 
  RCC2 %>%
  select(contains("mir499")) %>% 
  tidyr::pivot_longer(everything())
gghistogram(mir499_data, x = "value", facet.by = "name")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

mir499_data %>% 
  group_by(name) %>% 
  shapiro_test(value)

#not normal
#non-parametric * 2 paired groups --> Wilcoxon signed rank test

mir499_data %>% 
  wilcox_test(value ~ name, paired = TRUE,alternative = "greater")

mir499_data %>% 
  group_by(name) %>% 
  get_summary_stats(value)

#it increases after ttt
#ANSWER = REJECT   

3.5

gghistogram(RCC2, x = "DFFA", facet.by = "Grade")

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.

#not normal
#non-parametric * > 2 groups --> kruskal Wallis
RCC2 %>% 
  kruskal_test(DFFA ~ Grade)

#ANSWER= REJECT
#No need to perform post-hoc test