<table width = 900 border = 0><tr><td>
# <b>NMST539 | Lab Session 1</b>

## <b>Multivariate data visualisation using R</b>
### <font color = "#880000"><b>Summer Term 2021/2022 |  15/02/22 </b></font>
###### <a href = "Lab1_2022.Rmd">sourse file (UTF8 encoding)</a>

**************



Program R (GNU GPL public licence) can be downloaded from

<a href = "https://www.r-project.org">https://www.r-project.org</a>

There are suitable software distributions for OS Windows, Linux, as well as OS X (Macintosh}. 

Additional packages/libraries can be found at the repositary webpage <a href = "https://cran.r-project.org/mirrors.html">https://cran.r-project.org/mirrors.html</a>.
However, packages are usually programmed by ordinary R users and, therefore, some caution is always needed when using such libraries and the commands implemented in them. 
The functions/commands that we use in R can not serve as some software black boxes but, instead, the idea will be always placed on the underlying theoretical background (statistical 
theory and the probabilistic model) which should be always fully understood. 

User-friendly interfaces for the R software can be used in addition to simlify the work. One of the most common interface for R is the <a href = "https://www.rstudio.com">RStudio</a> 
(supports, again, OS Windows, Linux, and OS X/Macintosh).<br><br>

<b>Usefull references for statistical analyses in <i>R</i> (in Czech/in English)</b><br>
		
<ul>
<li>Bína, V., Komárek, A. a Komárková, L.: <i>Jak na jazyk R.</i>  (<a href = "NMFM301/Rmanual2.pdf">PDF súbor</a>)</li>
<li>Komárek, A.: <i>Základy práce s R.</i>  (<a href = "http://msekce.karlin.mff.cuni.cz/~komarek/vyuka/dataRko/Rmanual1.pdf">PDF súbor</a>)</li>
<li>Kulich, M.: <i>Velmi stručný úvod do R.</i> (<a href = "NMFM301/uvodr.pdf" >PDF súbor</a>)</li>
<li>De Vries, A. a Meys, J.: <i>R for Dummies.</i> (ISBN-13: 978-1119055808) </li>		
</ul><br>


### <font color = "#880000"><b>1. Grafical tools in R</b></font>

Three categories in particular: 

<ol>
<li><font color = "#880000"><b>Functions opening a new graphical device </b></font><br>
Each function from this category opens a new graphical device (window) which can be later used to add some 
addtitional elements into it.<br>

For instance: `plot()`, `dotchart()`, `hist()`, `barplot()`, `pie()`, `boxplot()`, `pairs()`, atď.;<br><br>
</li>

<li><font color = "#880000"><b>Functions that add elements into an active graphical window</b></font><br>
Functions from this category only work if there is already some graphical window openned (by using some function 
from the previous category). The functions can be used to add some additional plots/elements into an existing plot. <br>

For instance: `points()`, `lines()`, `abline()`, `text()`, `polygon()`, `rect()`, `box()`, `segments()`, `arrows()`, `legend()`, atď.;<br><br>
</li>

<li><font color = "#880000"><b>Functions meant for plots/graphs manipulation</b></font><br>
Functions from this category are primarily used for some manipulation with the existing/complete plot outputs -- such as 
saving or organizing them into some structure. <br>

For instance: `dev.new()`, `X11()`, `dev.off()`, `postscript()`, `pdf()`, `png()`, `par()`, `layout()`, atď.;<br><br>
</li>
</ol>

It is always good to use the R help which is always available for each R function. For each command, the corresponding help 
can be obtained by using the command itself given right after the question mark. For instnace, for the help for the command `plot()` 
one needs to type `?plot`. <br><br>


#### <font color = "#880000"><b>Scatterplots</b></font>

Simple, easily understandable 2-dimensional plots. However, it does not necessarily mean that two different covariates (two sets of information) can be only used to draw a scatterplot. 
Multivariate (three and more) dimensional scatterplots are, however, not recommended as human eyes are not well equiped for reading and understanding information in higher dimensions.

For illustration puposes we will use the dataset `mtcars` which is standardly available in the R software (uppon installation). 

```{r, eval = FALSE}
head(mtcars)
dim(mtcars)
summary(mtcars)
```

```{r, echo=FALSE}
attach(mtcars)
```

Simple scatterplot (using two covariates only) and the default version of the command `plot()`: 
```{r, echo=FALSE, fig.width = 10}
plot(mpg ~ cyl, data = mtcars)
```

The graph above contain only very limited information. For instance, it is not even clear from the plot itself what type of information is given in the figure. 
The golden rule says that <b>each figure, together with its caption (which should be necessary), must contain all relevant information needed to fully understand the plot</b>. 
Thus, lets try another plot -- containing the same information -- and much more: 

```{r, fig.width = 10}
x <- mtcars[order(mtcars$mpg),] # sort by mpg
x$cyl <- factor(x$cyl) # it must be a factor
x$color[x$cyl==4] <- "red"
x$color[x$cyl==6] <- "blue"
x$color[x$cyl==8] <- "darkgreen"  
dotchart(x$mpg,labels=row.names(x),cex=.7,groups= x$cyl,
  	main="Gas Milage for Car Models by cylinders",
   xlab="Consumption [Miles Per Gallon]", gcolor="black", color=x$color, pch = x$am)

c4 <- mean(mtcars$mpg[mtcars$cyl == 4])
c6 <- mean(mtcars$mpg[mtcars$cyl == 6])
c8 <- mean(mtcars$mpg[mtcars$cyl == 8])

lines(c(c8, c8), c(1, 14), col = "darkgreen", lwd = 2)
lines(c(c6, c6), c(17, 23), col = "blue", lwd = 2)
lines(c(c4, c4), c(26, 36), col = "red", lwd = 2)
```

Alternatíve with a matrix of scatterplots given by the command `pairs()`:

```{r, fig.width = 10}
pairs(mtcars[,c(1,3,4,6)], col=mtcars$cyl, pch = mtcars$am)
```


Or, instead, we can try some plots which already implicitly include some statistical methods 
(e.g., estimation of the densities or the corresponding regression lines between the given pairs) may be obtained by using the additional package `library(psych)` 
and the command  `pairs.panels()`:

```{r, message=F, warning=F, fig.width = 10}
library(psych)
pairs.panels(mtcars[,c("mpg", "disp", "hp", "wt")], 
             method = "pearson", # correlation method
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots
             ellipses = TRUE # show correlation ellipses
             )
```

There are of course many options how to create a good plot/graph using the R software. The main goal should be to offer 
an easily understandable graph with all relevant information being clearly contained in it. This also includes a proper description 
of the axes, indlucind the unites (for instance, meters for distances, liters for volumes, etc.).

An example of a plot with two different $y$ axes: 

```{r, fig.width = 10}
par(mar = c(5, 4, 4, 4) + 0.3)
x <- mtcars$wt
plot(mtcars$mpg ~ x, data = mtcars, pch = 21, bg = "lightblue", cex = 0.8, xlim = c(1, 5.5), xlab = "Total weight [1000 lbs]", ylab = "")
lines(lowess(mtcars$mpg ~ mtcars$wt), col = "black", lwd = 2)

mtext("Consumption [miles per gallon]" ,side=2,line=3, col = "black")
axis(2,col="black",col.axis="black")

par(new=TRUE)


plot(mtcars$hp ~ x, axes = FALSE, bty = "n", xlab = "", ylab = "", type = "p", col ="red", ylim = c(min(mtcars$hp), max(mtcars$hp)), lwd = 2, cex = 0.5, xlim = c(1, 5.5))
lines(lowess(mtcars$hp ~ x), col = "red", lwd = 2)
axis(4,col="red" ,col.axis= "red")

mtext("Horse power",side=4,line=3, col = "red")

for (i in 1:50){
  lines(c(1 + i * 4.5/50, 1 + i * 4.5/50), c(50, 340), col = "gray")
}
```



Sometimes, a very simple plot may just do the job which is needed... what can be seen in the following plot for instance? 

```{r, fig.width = 10}
library(MASS)
library("RColorBrewer")
colors <- brewer.pal(5, "Reds")

z <- kde2d(mtcars[,1], mtcars[,4], n=5)

plot(mtcars[,c(1,4)], xlab="Miles per Galon", ylab="Horse Power", pch=21, bg = "yellow")
contour(z, drawlabels=FALSE, nlevels=4, col=colors, add=TRUE)
```


The library `MASS` above is needed to calculate the kernel estimate for the two dimensional density (command `kde2d()`).  
In addition, the library `RColorBrewer` provides the colour scale for differed colors of the countours.  <br><br>

In practical applications the plots should always correspond with the question of interest and the main scientific hypothesis. 
If we are interested in a car consumption given the number of cylinders, mostly it will be completely irrelevant to give a plot 
of the car weights given the number of transmission gears.<br><br>


#### <font color = "#880000"><b>Boxplots</b></font>

Scatterplots are used for the visualisation of the raw data points. Instead, boxplots are used to visualise some useful descriptive characteristics -- sample
versions of some theoretical properties -- such as the sample median/mean, sample variance/tandard error, sample quartiles, etc.

The boxplot figure must be understood properly. It is fundamental for the interpretation. 

```{r, fig.width = 10}
boxplot(mpg ~ cyl, col = brewer.pal(3, "Blues"), xlab = "number of cylinders", ylab = "miles per galon", main= "")
```

Distinguishing for two different covariates (pay attention on the x axis): 

```{r, fig.width = 10}
boxplot(mpg ~ cyl + am, col = c(brewer.pal(3, "Greens"), brewer.pal(3, "Blues")), xlab = "No. of cylinders(.)Auto Transmission Indicator")
```

Or a sophisticated boxplot versions with the consumption given two covariates together with a 95% confidence interval for the true (unknown) median value (for `notched = TRUE`). 

```{r, fig.width = 10, warning = FALSE}
boxplot(mpg ~ cyl * am, col = c(brewer.pal(3, "Greens"), brewer.pal(3, "Blues")), xlab = "Number of Cylinders", ylab = "Miles per Galon | Automatic Transmission", main= "", notch = T)
```

Alternative  histograms  showing also the estimated density (using command `vioplot()` from the library `vioplot`). 
```{r, message=F, warning=F}
library(vioplot)
```

```{r, fig.width = 10}
x1 <- mtcars$mpg[mtcars$cyl==4]
x2 <- mtcars$mpg[mtcars$cyl==6]
x3 <- mtcars$mpg[mtcars$cyl==8]
vioplot(x1, x2, x3, names=c("4 cyl", "6 cyl", "8 cyl"), col="gold")
title("Violin Plots of Miles Per Gallon")
```

Similar boxplot can by also obtained by using the command  `pirateplot()` and the library  `yarrr`. For a detailed interpretation of the 
boxplot description, see <a href = "https://cran.r-project.org/web/packages/yarrr/vignettes/pirateplot.html">manuál</a>.

```{r, fig.width = 10, warning = F, message=F}
#install.packages("yarrr")
library("yarrr")
pirateplot(mpg ~ cyl, data = mtcars, xlab = "number of cylinders", ylab = "miles per galon", main= "", inf.method = "iqr")
```

THe same command, however, with the confidence intervals instead of the sample quartiles: 
(see again the <a href = "https://cran.r-project.org/web/packages/yarrr/vignettes/pirateplot.html">manuál</a>).

```{r, fig.width = 10, warning = F, message=F}
pirateplot(mpg ~ cyl, data = mtcars, xlab = "number of cylinders", ylab = "miles per galon", main= "", inf.method = "ci")
```



<br><br>

#### <font color = "#880000"><b>Parallel Coordinate Plots </b></font>

This type of plots may reveal some dependence form between various covariates which are given in the dataset. 
Again, the interpretation of the plot may be tricky in some sense but, sometimes, very useful information can 
be taken from the parallel coordinate plots. 

```{r, fig.width = 10, message = FALSE}
library(MASS)
attach(mtcars)
```

```{r, fig.width = 10, message = FALSE}
colorVector <- rep("brown", dim(mtcars)[1])
colorVector[am == 0] <- "brown1"
colorVector[cyl == 6] <- "blue4"
colorVector[cyl == 6 & am == 0] <- "blue1"
colorVector[cyl == 4] <- "darkgreen"
colorVector[cyl == 4 & am == 0] <- "green"
parcoord(mtcars, col = colorVector, lty = gear)
```

How would you interpret the plot above? What would be your conclussion? 


<br><br>

#### <font color = "#880000"><b>Library ggplot2 </b></font>

A very popular library for graphical outputs in R -- the R package `ggplot2`: 

```{r, fig.width = 10, message = FALSE, eval = FALSE}
install.packages("ggplot2")
```

```{r, fig.width = 10, message = FALSE}
library(ggplot2)
```

The main plotting function is  `qplot()` which, unlike the standard R command `plot()`, offers many more options how to create 
the final graph (see the help `?qplot()`).  

Let us slightly prepare the data from the `mtcars` object: 

```{r, fig.width = 10, message = FALSE}
mtcars$gear <- factor(mtcars$gear,levels=c(3,4,5), labels=c("3gears","4gears","5gears")) 
mtcars$am   <- factor(mtcars$am,levels=c(0,1), labels=c("Automatic","Manual")) 
mtcars$cyl  <- factor(mtcars$cyl,levels=c(4,6,8),  labels=c("4cyl","6cyl","8cyl"))
```

Density estimates for the car consumption given the number of geers (using `qplot()`):

```{r, fig.width = 10, message = FALSE}
qplot(mpg, data=mtcars, geom="density", fill=gear, alpha=I(.5), main="Distribution of Gas Milage", xlab="Miles Per Gallon",  ylab="Density")
```

The same function/command can be also used to get boxplots (using the additional parameter `geom = "..."`): 

```{r, fig.width = 10, message = FALSE}
qplot(gear, mpg, data=mtcars, geom=c("boxplot", "jitter"), fill=gear, main="Mileage by Gear Number", xlab="", ylab="Miles per Gallon")
```

Or even regression lines: 

```{r, fig.width = 10, message = FALSE, warning = FALSE}
qplot(wt, mpg, data=mtcars, geom=c("point", "smooth"), method="lm", formula=y~x, color=cyl, main="Regression of MPG on Weight", xlab="Weight", ylab="Miles per Gallon")
```   
   
<br><br>   
   
#### <font color = "#880000">Useful hint</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
The library `ggplot2` uses multiple layers to put together the final graphical output. The layers can be combined 
together using different commands, for instance, `geom_points()`, `geom_abline()`, `geom_area()` `geom_bar()`, `geom_contour()`, `geom_polygon()`. 
(see a detailed list at <a href = "http://sape.inf.usi.ch/quick-reference/ggplot2/geom">http://sape.inf.usi.ch/quick-reference/ggplot2/geom</a>). 

A simple illustration for multiple layers in one plot: 

```{r, fig.width = 10, message = FALSE, warning = FALSE}
dataSum <- plyr::ddply(mtcars, "gear", plyr::summarise, mean = mean(mpg), sd = sd(mpg))

ggplot() +
  geom_point(data = mtcars, aes(x = gear, y = mpg)) +
  geom_point(data = dataSum, aes(x = gear, y = mean), colour = 'red', size = 3) +
  geom_errorbar(data = dataSum, aes(x = gear, y = mean,  ymin = mean - sd, ymax = mean + sd), colour = 'red', width = 0.4) +
  ggtitle("Car consumption summary given the number of gears")
``` 

<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<hr size="3" width="100%" noshade style="color:#440000" align="left" /><br><br>   




<br>

In general, there are very many options in R to produce nice and (information) efficient plot/graph. All possibilities go
way beyond the scope of this lab session. Therefore, we can only refer a very nice manual with many practical examples 
which can be all found here 

<a href = "https://www.r-graph-gallery.com/">https://www.r-graph-gallery.com/</a> and also here <a href = "https://statisticsglobe.com/graphics-in-r">https://statisticsglobe.com/graphics-in-r</a>.

<br><br>



### <font color = "#880000"><b>2. Reporting results with Sweave and Knit</b></font>

Important part of each statistical analysis is the result reporting which also requires some substantial time amount. 
An appropriate form should be always used. There are some straightforward options in the R software. For instance, 
one can easily obtain a (very simply reconstructable) PDF report with all important analysis, relevant results, and the corresponding 
interpretation using the LaTeX and the R library `Sweave`. Alternatively, one can report the results in terms of some www document (HTML file) using the R library  `Knit`. 
The last one can be actually also used to create a PDF file. 

An illustrative guide can be found here: 

<a href = "https://support.rstudio.com/hc/en-us/articles/200552056-Using-Sweave-and-knitr"> https://support.rstudio.com/ </a>

The library/command `Sweave` enables cooperation between the R software and LaTeX. Analogously, the library/command `Knit` integrates the R software and the HTML code. 

The main principle relies on the  `tex` source file (for Sweave) or the  `html` source file (for Knit), which are saved as  `*.Rnw` (Sweave) or `*.Rmd` (Knit) types. 

In addition to a standard tex code (Sweave) or the html code (Knit) one can also use a (properly specified) R software code embeded in specific symbols. When applying the command  `Sweave("file_name.Rnw")` or `Knit(file_name.Rmd)` 
the R software reads the whole document firstly, the corresponding parts with the R code are evaluated and the corresponding output is plugged in instead. 
Consequently the R software uses LaTeX (Sweave) or HTML compiler (Knit) to produce the final PDF/HTML document/file.


The whole process simplifies when one uses the interface  `RStudio`

<a href = "https://www.rstudio.com"> https://www.rstudio.com </a>

where Sweave and Knit buttoms are available by default (they appear according to the type of the loaded source file -- `*.Rnw` (Sweave) or `*.Rmd` (Knit). 
A detailed guide can be found where: 

<a href = "http://gosset.wharton.upenn.edu/teaching/471/EPFL-Sweave-powerdot.pdf">http://gosset.wharton.upenn.edu/</a>

<br><br>

#### <font color = "#880000"><b>Sweave and Knit</b></font>

<ul>
<li>Illustrative example for the `Rnw` file:<br>
<a href = "sweave_test.Rnw">Download Rnw file</a> (Sweave | UTF 8 encoding)<br> 
The final PDF document is obtained when using the  `Compile PDF` buttom in the RStudio interface -- see the figure below.</li>

<li>Illustrative example for the `Rmd` file:<br>
<a href = "Lab1_2022.Rmd">Download Rnw file</a> (Knit | UTF 8 encoding)<br>
The final HTML document is obtained when using the `Knit HTML` buttom in the RStudio interface -- viď the figure below.</li><br><br>
</ul>


<center>
<img src = "rnw.png" width = "900">
<br><br><br>
<img src = "rmd.png" width = "900">
</center>
<br><br>

<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<hr size="3" width="100%" noshade style="color:#440000" align="left" /><br><br> 


## <font color = "#880000"><b>Individual task assignment</b></font>
### (Deadline: Second lab session / 22.02.2022)

A proper ellaboration of all partial tasks given at the end of each HTML Markdown (each week)
is an obligation for all students who intend to obtain the final credit at the end of the term. 

<ul>
<li>Use the guide on the web page of the  <a href = "https://labs.karlin.mff.cuni.cz/pocitacova-laborator-k10/">K10 lecture room</a> and create your own (simple) web page where you will 
publish all your solutions of the individual tasks 
assigned during the term (technically, it is sufficient to create one file called <i>index.html</i> which will be located in the directory <i>public_html</i> in your home directory). 

The file <i>index.html</i> can be edited using an arbitrary text editor or, instead, you can also edit the file directly in the R software (using the `Knit` library). 
Your web page link must be sent by emial to   <i>maciak [AT] karlin.mff.cuni.cz</i> and also to  <i>hlavka [AT] karlin.mff.cuni.cz</i>. </li><br>


<li>If you are interested, the page access can be also restricted by a password using <i>htaccess</i>.  In such case, you should also provide the corresponding password when sending the link according the the item above.</li>
</ul>

Use the information provided on <a href = "https://www2.karlin.mff.cuni.cz/~hlavka/teac.html">the lecturer's website </a> to install the R library `SMSdata` with various 
datasets inclded in this packages. From all the datasets available in the library `SMSdata`


```{r, eval = F}
data(package = "SMSdata") 
```
 
 choose one specific dataset  (free choice, there are 21 of them in total). The choosen dataset can be loaded into the R environment by using the command
`data(dataset_name)`, for example, 

```{r, eval = F}
data(plasma) 
```

<ul>
<li>Choose one dataset from the `SMSdata` library and load into R;</li>

<li>Use some proper graphical tools to visualize (primarily the multivariate character of the) data;</li> 

<li>Upload at least two figures on your web site. Add some reasonable interpretation/comments and discuss (at least briefly)
the multivariate nature of the data using the created plots; </li>

<li>The web site link with the figures and comments should be sent by email to <i>maciak [AT] karlin.mff.cuni.cz</i> and <i>hlavka [AT] karlin.mff.cuni.cz</i> on
<font color = "red">Tuesday, 21.02.2022, at 12:00 at latest</font>.</li>
</ul>

</td></tr></table>
<br><br><br>