# 8 Bivarijantno modeliranje: dvije kontinuirane varijable, 1. dio


```{r}
getwd()
```

```{r}
setwd("rossiter/scripts/")
```

```{r}
mhw = read.csv("../mhw.csv")
```

Pogledajmo strukturu data framea `mhw`.

```{r}
str(mhw)
```

Grafički prikažimo točke sa kordinatama `grain` i `straw`.

```{r fig.width=7, fig.height=6, tidy=FALSE}
plot(mhw$straw ~ mhw$grain,
     xlab = "Prinos zrna",
     ylab = "Prinos slame",
     pch = 20, col="darkblue", cex=1.2,
     sub="Median je prikazan pomocu crvenih isprekidanih linija")
title (main="Slama vs. zrnevlje, Mercet-Hall eksperiment")
grid()
abline(v=median(mhw$grain), lty=2, col="red")
abline(h=median(mhw$straw), lty=2, col="red")
```

Nacrtajmo "istu" sliku pomoću paketa *ggplot2*.

```{r fig.width=7, fig.height=6, tidy=FALSE}
ggplot(mhw, aes(grain, straw)) +   
  geom_point(col="darkblue", size=3) +   
  ggtitle("Slama vs. zrnevlje, Mercet-Hallov eksperiment") +   
  geom_hline(aes(yintercept=median(straw)), col="red", lty=2) +   
  geom_vline(aes(xintercept=median(grain)), col="red", lty=2) +   
  xlab("Prinos zrna\nCrveno je naznacen polozaj mediana") + ylab("Prinos slame") +   
  theme_bw()+   
  theme(plot.title = element_text(face="bold", colour="#990000", size=20))
```

## Korelacija zrnevlja i slame

Dva su najznačajnija mjerila korelacije

- Parametarski: Pearsonov koeficijent korelacije $r$ i
- Neparametarski: Spearmanov koeficijent $\rho$.

Ovdje ćemo upotrebiti Pearsonov koeficijent. Ako je dvodimenzionala
razdioba normalna tada taj koeficijent ima neka dobra svojstva.
To je procjenitelj korelacije populacije po metodi maksimalne
vjerodostojnosti, zatim je asimptotski nepristrasan i efikasan procjenitelj.

Stoga ćemo naša mjerenja usporediti sa uzorkom dvodimenzionalne
normalne razdiobe sa istim parametrima kao što su 
parametri empirijske razdiobe.

### Generiranje dvodimoenzionalne normalne razdiobe i usporedba sa mjerenim podacima

**Zadatak**: Vizualiziraj dvodimenzionalnu normalnu razdiobu sa parametrima  iz prinosa zrna i slame; vizualno usporedi uzorak sa podacima iz `mhw`.

Za dobivanje podataka dvodimenzionalne korelirane normalne razdiobe
koristiti ćemo funkciju `mvrnorm` iz biblioteke `MASS`. Za to nam 
trebaju podaci o srednjim vrijednostima i *variance-covariance* matrica.

Vrijednosti srednjih vrijednosti ćemo dobiti pomoću funkcije `colMeans()`.

```{r}
mu = colMeans(mhw[, c("grain", "straw")])
mu
```

*Varijace-covariance* matricu ćemo dobiti pomoću vektorizirane `var()` 
funkcije

```{r}
sigma = var(mhw[, c("grain", "straw")])
sigma
```

Sada mižemo generirati uzorak iz dvodimenzionalne 
korelirane normalne razdiobe.

```{r}
require(MASS)
sim.sample= mvrnorm(length(mhw$grain), mu = mu, Sigma = sigma)
head(sim.sample)
```

Usporedimo osnovne podatke genriranih i mjerenih veličina.
```{r}
summary(sim.sample)
```

```{r}
summary(mhw[, c("grain", "straw")])
```

####  Grafički usporedimo rezultate mjeranja i simulirane podatke.

##### U prvoj usporedbi ćemo koristiti histograem.

```{r, tidy=FALSE}
par(mfrow = c(2, 2))
grain.lim = c(min(sim.sample[, "grain"], mhw$grain ), 
              max(sim.sample[, "grain"], mhw$grain))
straw.lim = c(min(sim.sample[, "straw"], mhw$straw), 
              max(sim.sample[, "straw"], mhw$straw))
hist(mhw$grain, xlim = grain.lim, main = "Grain (actual)",
     col = "lightyellow", breaks = seq(grain.lim[1], grain.lim[2],
                                      length = 17))
hist(sim.sample[, "grain"], xlim = grain.lim, main = "Grain (simulated)",
     col = "antiquewhite", breaks = seq(grain.lim[1], grain.lim[2],
                                      length = 17))
hist(mhw$straw, xlim = straw.lim, main = "Straw (actual)",
     col = "lightblue", breaks = seq(straw.lim[1], straw.lim[2],
                                     length = 17))
hist(sim.sample[, "straw"], xlim = straw.lim, main = "Straw (simulated)",
     col = "lavender", breaks = seq(straw.lim[1], straw.lim[2],
                                    length = 17))
par(mfrow = c(1, 1))
```

Nacrtajmo "istu" sliku pomoću bibilioteke *ggplot2*.

Za crtanje više figura na jednoj slici treba nam funkcija `multiplot()`

```{r, tidy=FALSE}
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
```


```{r, tidy=FALSE}
ssdf = data.frame(sim.sample)

grain.lim = c(min(sim.sample[, "grain"], mhw$grain ), 
              max(sim.sample[, "grain"], mhw$grain))
straw.lim = c(min(sim.sample[, "straw"], mhw$straw), 
              max(sim.sample[, "straw"], mhw$straw))

binwidth.grain = (grain.lim[2] - grain.lim[1])/16
binwidth.straw = (straw.lim[2] - straw.lim[1])/16
p1 = ggplot(mhw, aes(grain)) + 
  geom_histogram(binwidth=binwidth.grain, col="black", fill="lightyellow") +
  coord_cartesian(xlim = grain.lim) +
  ggtitle("Grain (actual)") + theme_bw() + 
  theme(panel.grid.major=element_line(size = 0.3, color="black",
                                      linetype = "dotted"))
p2 = ggplot(ssdf, aes(grain)) +
  geom_histogram(binwidth=binwidth.grain, col="black", fill="antiquewhite") +
  coord_cartesian(xlim = grain.lim) +
  ggtitle("Grain (simulated)") + theme_bw()+ 
  theme(panel.grid.major=element_line(size = 0.3, color="black",
                                      linetype = "dotted"))
p3 = ggplot(mhw, aes(straw)) + 
  geom_histogram(binwidth=binwidth.straw, col="black", fill="lightblue") +
  coord_cartesian(xlim = straw.lim) +
  ggtitle("Straw (actual)") + theme_bw()+ 
  theme(panel.grid.major=element_line(size = 0.3, color="black",
                                      linetype = "dotted"))
p4 = ggplot(ssdf, aes(straw)) +
  geom_histogram(binwidth=binwidth.straw, col="black", fill="lavender") +
  coord_cartesian(xlim = straw.lim) +
  ggtitle("Straw (simulated)") + theme_bw() + 
  theme(panel.grid.major=element_line(size = 0.3, color="black",
                                      linetype = "dotted"))

multiplot(p1, p2, p3, p4, cols=2)
```

# Formule u LaTeX-u

$$
\phi(x) = \int_0^x \frac{t\,dt}{1+\sin^3(t)}
$$


Tables
======
 
```{r table1, comment=NA, results='asis'}
library(xtable)
print(xtable(head(mhw, 6)), type = "html", include.rownames = F)
```

- Use comment=NA in the knitr chunk options to prevent lines of output beginning with ##
- Use results=’asis’ in the knitr chunk options so as the HTML is rendered, not printed
- Use type = “html” as an argument to print xtable() for HTML tables



##### Za drugu usporedbu koristiti ćemo grafički prikaz točaka

```{r fig.width=10, fig.height=6, tidy=FALSE}
par(mfrow = c(1, 2))

plot(sim.sample, main = "Slama vs. zrnevlje (simulirano)", 
     xlab = "Zrnevlje", ylab="Slama",
     xlim = grain.lim, ylim = straw.lim,
     pch = 20, col= "blue" )
abline(v = median(sim.sample[, 1]), lty = 2, col = "red")
abline(h = median(sim.sample[, 2]), lty = 2, col = "red")
plot(mhw$grain, mhw$straw, main = " Slama vs. zrno (izmjereno)",
    xlab = "Zrnevlje", ylab="Slama",
    xlim = grain.lim, ylim = straw.lim,
    pch = 20, col= "black" )
abline(v = median(mhw$grain), lty = 2, col = "red")
abline(h = median(mhw$straw), lty = 2, col = "red") 

par(mfrow = c(1,1))
```

Analognu sliku nacrtajmo pomoću biblioteke *ggplot2*.

```{r fig.width=10, fig.height=6, tidy=FALSE}
p1 = ggplot(ssdf, aes(grain, straw)) +
  geom_point(color = "blue") +
  labs(title="Slama vs zrno (simulirano)", x = "Zrno", y = "Slama") +
  geom_hline(aes(yintercept=median(straw)), col="red", lty=2) +
  geom_vline(aes(xintercept=median(grain)), col="red", lty=2) +
  theme_bw() +
  theme(panel.grid.major=element_line(size = 0.3, color="black",
                                      linetype = "dotted"),
        plot.title = element_text(face="bold", colour="#990000", size=15),
        axis.title.x = element_text(face="bold", size=12),
        axis.title.y = element_text(face="bold", size=12))

p2 = ggplot(mhw, aes(grain, straw)) +
  geom_point(color = "black") +
  labs(title="Slama vs zrno (izmjereno)", x = "Zrno", y = "Slama") +
  geom_hline(aes(yintercept=median(straw)), col="red", lty=2) +
  geom_vline(aes(xintercept=median(grain)), col="red", lty=2) +
  theme_bw() +
  theme(panel.grid.major=element_line(size = 0.3, color="black",
                                      linetype = "dotted"),
        plot.title = element_text(face="bold", colour="#990000", size=15),
        axis.title.x = element_text(face="bold", size=12),
        axis.title.y = element_text(face="bold", size=12))

multiplot(p1, p2, cols=2)
```

Budući da smo gotovi sa vizualnom usporedbom izbrisati cemo nepotrebne objekte.

```{r}
rm (sim.sample, grain.lim, straw.lim, binwidth.grain, binwidth.straw)
```

### Koeficijent korelacija za mjerene i simulirane podatke.


Kako se iz predhodnog vidi empirijska dvodimenzionalna razdioba
zrna i slame je približno normalna pa će
Pearsonov koeficijent korelacije imati mnoga dobra svojstva.

**Zadatak** Izračunajmo koeficijent korelacije između zrnevlja i slame i njegov 95% interval povjerenja.

```{r}
cor (mhw$grain, mhw$straw)
```

Usporedimo dobiveni koeficijent korelacija sa rezultatom simuliranih podataka.

```{r}
cor (ssdf$grain, ssdf$straw)
```

Vidimo da su dobiveni rezultati slični.

Odredimo interval povjerenja korelacije.

```{r}
rez = cor.test(mhw$grain, mhw$straw)
rez
```

Ako nas zanima samo interval povjerenja a ne cjelokupni output funkcije `cor.test()`
izvedemo naredbu
```{r}
rez$conf.int
```