---
title: "Introduction to Biostatistics <br> The 1st Lecture"
author: "Hiroyoshi IWATA  hiroiwata@g.ecc.u-tokyo.ac.jp"
date: "2022/4/8"
output:
  html_document: default
---


## Packages

Packages required for executing this Rmd code.

```{r, message = F}
require(KernSmooth)
require(maps)
require(mapdata)
require(plotly)
```

## Introduction

  Recently, in the fields of agriculture and life sciences, various types of data are being collected and accumulated in large quantities. It is necessary to analyze the data using methods appropriate to the purpose of the research and the nature of the data, in order to ensure that we do not miss the “undiscovered knowledge” harbored in the data.

  There are various methods for statistical analysis. In order to understand the features of the methods and the principles of analysis and also to interpret the analysis results properly, you should learn them accordingly. Also, in order to make your learning more effective, it is essential to have the experience of analyzing actual data. It is often the case that you can clearly understand what you have learned in lectures and textbooks when you analyze your data by yourself.

  This course is intended to provide you with a "first step" to analyze their own data and improve their statistical analysis skills. Specifically, we will focus on practical data analysis methods using R, focusing on some of the statistical methods that may be required in your future research. The goal of this lecture is to acquire the skills to use it for data analysis of general-purpose statistical analysis methods, e.g., regression analysis, analysis of variance, and principal component analysis. Furthermore, the goal is to build a foothold to perform more advanced data analysis. Although it is a short course with four lectures, I will provide you this course so that you can be interested in the joy and skills of statistical analysis.

## R
  R is free software for statistical analysis. (To put it a little correctly, R is the name of a computer language. R installed on a PC as software is an "environment" for using R language). R has many functions, and its usages range from statistical analysis to pre-processing of data, overview of data, and creation of graphs for papers. In addition, various analysis can be easily performed by installing an extension program distributed as a “package”. Newly developed statistical methods will be available quickly in R. Thus, the skills for using R have become useful to researchers in agronomics and life sciences.

In addition, for learning R, a large number of reference books are available. Introductory books I recommend are:

- Peter Dalgaard, Introductory Statistics with R (Statistics and Computing) Second Edition, Springer, 2008, ISBN: 978-0387790534

- Brian Everitt, Torsten Hothorn, An Introduction to Applied Multivariate Analysis with R (Use R!), Springer, 2011, ISBN: 978-1441996497

## Simple calculation using R
  In analysis using R, the analysis is basically progressed interactively while sequentially inputting commands (However, when you actually perform analysis, it is useful to prepare a series of commands, as an R script, prior to the analysis. Remember that it is more convenient to execute the R script, because it allows us to do partial corrections and to review the history of analysis).

  Let's start with getting used to R while doing simple calculations with command input.

The easiest way to use R is to enter a simple arithmetic expression and get the answer. For example,
```{r}
3 + 5 * 3
```

If you want to perform the next calculation based on the obtained result, assign the resulted value to some variable as follows.
```{r}
x <- 1 + 2
x
```

The assigned value can be used for another calculation through the variable name.
```{r}
x + 5 * x
```

Various calculations can be performed using functions.
```{r}
abs(x)
sin(x)
atan(x)
log(x)
log10(x)
```

Let's perform a bit more complicated calculation. The probability density function of the normal distribution of mean $\mu$ and the variance $\sigma^2$ is
$$f(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$
Calculate this with R.
```{r}
mu <- 3
s2 <- 2
x <- 5
1 / sqrt(2 * pi * s2) * exp(- (x - mu)^2 / (2 * s2)) 
```

If you use the “dnorm” function to calculate the probability density of the normal distribution for confirming your calculation, you can get the same result.
```{r}
dnorm(x, mu, sqrt(s2))
```

```{r, echo = F, fig.cap = "Fig. 1. Normal distribution with mean 3 and variance 2"}
x <- seq(-8, 12, 0.1)
plot(x, dnorm(x, mu, sqrt(s2)), type = "l")
```

### Quiz1
  Now, let's solve a practice question here. Practice questions will be presented in the lecture.

Go to https://www.menti.com/d9p29idxtw , then, register your nickname and wait for the quiz to start.

## Calculation using vector or matrix
　A great advantage of R is that we can perform vector and matrix operations easily. Let's calculate some summary statistics using vector and matrix operations.

For example, we can easily create a vector of six elements as follows: this is the data which measured the grain length of six varieties/lines of rice in mm unit (the source of the data will be described later).
```{r}
length <- c(8.1, 7.7, 8.2, 9.7, 7.1, 7.3)    # mm scale
length
```

Input grain widths of the same varieties/lines, and calculate the length-width ratio.
```{r}
width <- c(3.7, 3.0, 2.9, 2.4, 3.3, 2.5)
ratio <- length / width
ratio
```

First, let's calculate the average of the length-width ratio of grains. The estimate of the population mean is
$$\sum_{i = 1}^{n} x_i / n$$
where $x_i$ is the value of the $i$th sample and $n$ is the number of samples.
```{r}
sum(ratio)
length(ratio)
sum(ratio) / length(ratio)
```

The mean can be calculated using the “mean” function.
```{r}
mean(ratio)
```

Next, let's calculate the variance. An estimate of the population variance is
$$\sum_{i = 1}^{n} (x_i - \overline{x})^2 / (n-1)$$
where $\overline{x}$ is the average calculated earlier.
```{r}
xbar <- mean(ratio)
(ratio - xbar)^2
sum((ratio - xbar)^2)
sum((ratio - xbar)^2) / (length(ratio) - 1)
```

The variance can be calculated using the “var” function.
```{r}
var(ratio)
```

Next, let's calculate the covariance. An estimate of the covariance between bivariate x and y is
$$\sum_{i = 1}^{n} (x_i - \overline{x})(y_i - \overline{y}) / (n-1)$$
where $\overline{x}$ and $\overline{y}$ represent the mean of each variable.

```{r}
xbar <- mean(length)
ybar <- mean(width)
sum((length - xbar) * (width - ybar)) / (length(length) - 1)
```

The variance can be calculated using the “var” function.
```{r}
cov(length, width)
```

Let's calculate Pearson product moment correlation coefficient (hereinafter referred to just as correlation coefficient). The correlation coefficient is 
$$\frac{\sum_{i = 1}^{n} (x_i - \overline{x})(y_i - \overline{y})} {\sqrt{\sum_{i = 1}^{n}(x_i - \overline{x})^2}\sqrt{\sum_{i = 1}^{n}(y_i - \overline{y})^2}}$$

```{r}
s12 <- sum((length - xbar) * (width - ybar))
s1 <- sum((length - xbar)^2)
s2 <- sum((width - ybar)^2)
s12 / (sqrt(s1) * sqrt(s2))
```

As you can see in the equation, the correlation coefficient is the covariance divided by the standard deviation of both variables. Let's actually calculate it and check the result.
```{r}
cov(length, width) / (sd(length) * sd(width))
```

The correlation coefficient is standardized by dividing it by the standard deviation of both variables. Unlike the covariance, we can understand the relationship between variables without being influenced by the scale of the measurement value. Thus, it is suitable for comparing the strength of relationships between variables measured at different scales (such as weight and length).

Note that we can also calculate the correlation coefficient using the “cor” function.
```{r}
cor(length, width)
```

Now let's calculate variance and covariance using matrix calculations. First, combine length and width to create a 6 by 2 matrix.
```{r}
x <- cbind(length, width)
x
```

Then we use the function apply to find the average of each column.
```{r}
m <- apply(x, 2, mean)
m
```

Subtract the column average from each column.
```{r}
z <- sweep(x, 2, m)
z
```

After that, we can calculate variance and covariance (variance-covariance matrix) by using the product of the matrix.
```{r}
t(z) %*% z / (nrow(z) - 1)
```
The diagonal components are variances, and the nondiagonal components are covariances.

The variance-covariance matrix can be calculated with the cov function.
```{r}
cov(x)
```

### Quiz2
  Now, let's solve a practice question here. Practice questions will be presented in the lecture.

If you close the page of the quiz, go to https://www.menti.com/d9p29idxtw

## Import and analyze external data
　If you use R for your own research, I think that most of the time you read and analyze data organized with spreadsheet software. Here, I will explain the procedure for reading data prepared by other software into R and analyzing it. Here, we will use the data which was used in genome-wide association studies of rice genetic resources (Zhao et al. 2011; Nature Communications 2: 467). The data can be downloaded from Rice Diversity ( http://www.ricediversity.org/data/ ). 

The function “read.csv” is used to read the file saved in csv format.
```{r}
pheno <- read.csv("RiceDiversityPheno.csv", stringsAsFactors = T)
```

To check the size of imported data or a part of the data as follows.
```{r}
dim(pheno)
head(pheno[,1:4])
```

The data have a separate file that describes the origin of each genetic resource. Here, we load the file and combine it with “pheno” data. First, read the file.
```{r}
line <- read.csv("RiceDiversityLine.csv", stringsAsFactors = T)
head(line)
```

Since NSFTV.ID in the “line” data and NSFTVID in the “pheno” data correspond to each other, the two data are combined based on the information of these columns.
```{r}
data <- merge(line, pheno, by.x = "NSFTV.ID", by.y = "NSFTVID")
head(data[,1:14])
```

### Quiz3
  Now, let's solve a practice question here. Practice questions will be presented in the lecture.

If you close the page of the quiz, go to https://www.menti.com/d9p29idxtw

## Analysis of the data
　We often want to analyze a large number of variables to look at their distributions and relationships among them. Measurement data, howover, often include missing entries with some experimental reasons. Here, we will analyze the data which we load from the csv file.

Let's calculate the length-width ratio of grains and its average in the same way as mentioned above.
```{r}
ratio <- data$Seed.length / data$Seed.width
mean(ratio)
```
We cannot calculate the average, and only get the value of “NA”. Why is that?

This is because the ratio contains missing values (represented as NA in R).
```{r}
ratio[1:14]
```

In such cases, specify the option na.rm in the calculation.
```{r}
mean(ratio, na.rm = T)
```

Find the mean for all variables in the data with the “sapply” function.
```{r}
sapply(data[, 1:14], mean, na.rm = T)
```
A warning message will be displayed for data that is not numeric data, and the calculated result will be NA.

Using the following command, we can calculate not only averages but also quartiles, minimum values and maximum values for numeric data, and can count up samples belonging to each class for factor data. 
```{r}
summary(data[,1:14])
```

Let's calculate the correlation coefficient for grain length and width.
```{r}
cor(data$Seed.length, data$Seed.width)
```
The result is NA. This is due to missing data as before.

Specify the option for dealing with the missing value and try to calculate again.
```{r}
cor(data$Seed.length, data$Seed.width, use = "pair")
```
This time, the correlation is calculated successfully.

### Quiz4
  Now, let's solve a practice question here. Practice questions will be presented in the lecture.

If you close the page of the quiz, go to https://www.menti.com/d9p29idxtw

## Data visualization
　It is very important to look at the data from different angles before actually performing statistical analysis. For example, statistics such as mean and variance mentioned above are statistics for summarization, and even with variables with similar mean and variance, the distribution of observed values may differ greatly. Therefore, looking at the data carefully is important to understand the characteristics of the data. Data visualization is also necessary when we prepare the results of the analysis for the publication of a paper. Here we will learn various data visualization techniques.

First, let's make it possible to directly call variables in the data before explaining the visualization methods.
```{r}
attach(data)
```
With the “attach” command, for example, we can now specify Plant.height without data\$~. Otherise we have to specify the variable data\$Plant.height, 

Let's draw a histogram first.
```{r}
hist(Plant.height)
```

Let's draw a stem-and-leaf plot.
```{r}
stem(Plant.height)
```
This is not a graph. The result is shown with text.

Let's draw a box plot.
```{r}
boxplot(Plant.height)
```

Next, we will draw a histogram of blast resistance (Blast.resistance).
```{r}
hist(Blast.resistance)
```
It looks like the distribution is visualized well, but there is a “pitfall”.

The resistance to blast disease is expressed by the level of resistance with a score of 9 (0-9). Therefore, let's summarize how many accessions are included in each of 9 levels.
```{r}
t <- table(Blast.resistance)
t
```
You can see that the histogram we drew did not represent the whole class well.

You can draw a bar plot from the data summarized using the “table” function as described above.
```{r}
plot(t)
plot(t, xlab = "Blast resistance scores", ylab = "Frequency")
```

A bar chart can also be drawn using the “barplot” function. However, it looks a bit different from the bar chart drew above. 
```{r}
barplot(t)
```

Drawing a pie chart allows you to illustrate the percentage of each score.
```{r}
pie(t)
pie(t, main = "Blast resistance")
```

From here, let's look at the relationship between two variables.
```{r}
plot(Plant.height, Panicle.length)
```

Fit a straight line to the data by regression analysis.
```{r}
plot(Plant.height, Panicle.length)
abline(lm(Panicle.length ~ Plant.height))
```

Overlap the rug (textile) plot. This is useful for visualizing distribution density.
```{r}
plot(Plant.height, Panicle.length)
abline(lm(Panicle.length ~ Plant.height))
rug(Plant.height, side = 1)
rug(Panicle.length, side = 2)
```

Let's illustrate the relationship between the two variables using kernel smoothing.
```{r}
library("KernSmooth")
x <- data.frame(Plant.height, Panicle.length)
x <- na.omit(x)
d <- bkde2D(x, bandwidth = 4)
plot(x)
contour(d$x1, d$x2, d$fhat, add = T)
```
The contours represent the density of the points smoothed by the kernel.

Let's display this smoothed density in three dimensions.
```{r}
persp(d$x1, d$x2, d$fhat, xlab = "Plant.height", ylab = "Panicle.length", zlab = "density", theta = -30, phi = 30)
```

The data of Zhao et al. (2011) read from the file contain not only trait data but also data of the genetic background of genetic resources. Let's visualize both data together and investigate what kind of relationship between genetic background and trait.

SThe variable Sub.population represents differences in the genetic background of genetic resources. This is estimated using Structure analysis (Pritchard et al. 2000, Genetics 155: 945). Now, let's visualize and see what kind of relationship between genetic background and plant height and ear length.
```{r}
pop.id <- as.numeric(Sub.population)
plot(Plant.height, Panicle.length, col = pop.id)
levels(Sub.population)
legend("bottomright", levels(Sub.population), col = 1:nlevels(Sub.population), pch = 1)
```

Let's show in boxplots how there are differences in values due to differences in genetic background.
```{r}
boxplot(Plant.height ~ Sub.population)
boxplot(Plant.height ~ Sub.population, border = 1:nlevels(Sub.population))
```

Let’s see the relationships among plant height, panicle length, and flag leaf length with a bubble plot. Here, the size of the bubble represents the flag leaf length.
```{r}
symbols(Plant.height, Panicle.length, circles = Flag.leaf.length, inches = 0.1, fg = pop.id)
```

Let's draw a scatterplot of the relationship between the three variables for all possible combinations.
```{r}
x <- data.frame(Plant.height, Panicle.length, Flag.leaf.length)
pairs(x, col = pop.id)
```

Let's make it a bit more complex scatterplot through adding regression lines.
```{r}
pairs(x, panel = function(x, y, ...) {
	points(x, y, ...)
	abline(lm(y ~ x), col = "gray")
}, col = pop.id)
```

The data also have the latitude and longitude of the place where each genetic resource originates. Let's map and confirm the origin of each genetic resource on the world map.
```{r}
library(maps)
library(mapdata)
map('worldHires')
points(Longitude, Latitude, col = pop.id)
legend("bottomleft", levels(Sub.population), col = 1:nlevels(Sub.population), pch = 1)
```

The above command only draws a much smaller number of points than the number of genetic resources. This is because genetic resources from the same area overlap each other. In order to prevent the overlapping, move overlapped points a little with the function “jitter”.
```{r}
map('worldHires')
points(jitter(Longitude, 200), Latitude, col = pop.id)
legend("bottomleft", levels(Sub.population), col = 1:nlevels(Sub.population), pch = 1)
```

### Quiz5
  Now, let's solve a practice question here. Practice questions will be presented in the lecture.

If you close the page of the quiz, go to https://www.menti.com/d9p29idxtw

## Output to file of diagram
　It is useful to be able to output a graph to a PDF file in order to use the graph in a paper or presentation. Here, the method to do that will be explained briefly.

Let's output the figure we drew earlier to a file called map.pdf.
```{r}
pdf("map.pdf")
map('worldHires')
points(jitter(Longitude, 200), Latitude, col = pop.id)
legend("bottomleft", levels(Sub.population), col = 1:nlevels(Sub.population), pch = 1)
dev.off()
```
When the above command is executed, a file called map.pdf is output to the working directory of R.

Function “pdf” can specify the size of the output graph. When a larger size is required as is this figure, it is better to specify the size of the graph output explicitly.
```{r}
pdf("map_large.pdf", width = 20, height = 10)
map('worldHires')
points(jitter(Longitude, 200), Latitude, col = pop.id)
legend("bottomleft", levels(Sub.population), col = 1:nlevels(Sub.population), pch = 1)
dev.off()
```

In addition, if multiple figures are repeatedly output to the same pdf file, they will be saved in a pdf file with multiple pages. If you want to output similar and a large number of figures, it may be convenient to combine them into a single pdf file.

## Draw an interactive figure
You can use the package “plotly” to draw interactive diagrams. Here, let's draw the figure drawn earlier using plotly function plot_ly.

Try to display the density of data in three dimensions.
```{r, eval = FALSE}
require(plotly)
plot_ly(data = data.frame(d), x = d$x1, y = d$x2, z = d$fhat) %>%
add_surface()
```

Finally, we look at the relationship between the three variables.
```{r, eval = FALSE}
df <- data.frame(Sub.population, Plant.height, Panicle.length, Flag.leaf.length)
df <- na.omit(df)
plot_ly(data = df, x = ~Plant.height, y = ~Panicle.length, z = ~Flag.leaf.length, color = ~Sub.population, type = "scatter3d", mode = "markers")
```

### Quiz 6
  Now, let's solve a practice question here. Practice questions will be presented in the lecture.

If you close the page of the quiz, go to https://www.menti.com/d9p29idxtw

# Report assignment
Use at least four types of graphs learned in the lecture to illustrate the relationship between traits and/or the relationship between traits and genetic background. Describe the relationships that can be read from the graphs (figures) that you drew.

Submission method:

* Create a report as a pdf file and submit it to ITC-LMS.
* When you cannot submit your report to ITC-LMS with some issues, send the report to report@iu.a.u-tokyo.ac.jp
* Make sure to write the affiliation, student number, and name at the beginning of the report.
* The deadline for submission is April 22nd.

