---
title: "Multilivel Models -- Radon"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Reading in the data

```{r message=F}
library(ggplot2)
library(arm)
```

```{r}
setwd("/media/guerzhoy/Windows/sml310/W10/303/lec9/")
srrs2 <- read.table ("srrs2.dat", header=T, sep=",")
mn <- srrs2[srrs2$state=="MN",]
mn$radon <- mn$activity
mn$log.radon <- log (ifelse (mn$radon==0, .1, mn$radon))
```


## Complete pooling

```{r}
mean(mn$log.radon)
sd(mn$log.radon)
```

Let's compute the average log-radon for all counties



```{r}
fit <- lm(log.radon~county, data=mn)
s <- summary(fit)
coefs <- as.data.frame(s$coefficients)
```

Add a column with County nmaes


```{r}
cleanup_name <- function(s){
  return(gsub(" ", "", gsub("county", "", s)))
}


coefs$county <- as.character(lapply(rownames(s$coefficients), cleanup_name))
```


```{r}

ggplot(coefs, aes(x = county, y = Estimate)) +  
  geom_errorbar(aes(ymin=Estimate-`Std. Error`, ymax=Estimate+`Std. Error`))

subset_coefs <- coefs[c(2, 10, 22, 30, 36, 40, 50, 55, 60, 65),]
ggplot(subset_coefs, aes(x = county, y = Estimate)) +  
  geom_errorbar(aes(ymin=Estimate-`Std. Error`, ymax=Estimate+`Std. Error`))

```  
  
How many datapoints for Lac Qui Parle?

```{r}
mn$county_clean <- as.character(lapply(mn$county, cleanup_name))
mn[mn$county_clean=="LACQUIPARLE", c("zip", "log.radon")]
```

Just two datapoints. 

What about Anoka?

```{r}
mn[mn$county_clean=="ANOKA", c("zip", "log.radon")]
```

Lots more. In general, larger error bars means smaller sample size. We should trust that estimate less.


## Multilvel Model -- fake data

Let's generate county data

```{r}
mu_a <- 1.1
sigma_a <- .3
sigma_y <- 0.8

n_counties <- 10
county_alphas <- rnorm(n_counties, mu_a, sigma_a)
county_alphas
```

Now, 10 measurements for county 2 would look like this:

```{r}
rnorm(10, county_alphas[2], sigma_y)
```

Another 15 measurements for county 8 would look like this:

```{r}
rnorm(15, county_alphas[8], sigma_y)
```




## Partial Pooling

Let's now fit the Partial Pooling model

```{r}
library(lme4)
summary(lmer(log.radon~(1|county), data=mn))
```

Let's get the individual intercepts

```{r}
coef(lmer(log.radon~(1|county), data=mn))
```

Here, we are saying that we want to fit an individual intercept for each county, and each individual intercept for each county is modelled: that is, we are assuming that $\alpha_j~N(\mu_\alpha, \sigma_\alpha)$.



Let's look at Lac Qui Parle County specifically

```{r}
#Partial pooling
coef(lmer(log.radon~(1|county), data=mn))$county["LAC QUI PARLE       ",]

#No-pooling
as.data.frame(lm(log.radon~county, data=mn)$coefficients)["countyLAC QUI PARLE       ",]

#Complete pooling  
lm(log.radon~1, data=mn)
```

As you can see, partial pooling helps move the estimate towards the mean, by using information about the other counties.


## Adding predictors

Let's add the floor predictor

```{r}
lmer(log.radon~(1|county) + floor, data=mn) 
```

Let's fit a different intercept for each county as well.

```{r}
lmer(log.radon~floor+(floor|county) , data=mn) 
```

Let's get the modelled/random effects
```{r}
fit<-lmer(log.radon~floor+(floor|county) , data=mn) 
ranef(fit)
```


Now, let's do prediction. Suppose we want to make a prediction for a new measurement in Lac Qui Parle county.

```{r message=F}
fit<-lmer(log.radon~floor+(floor|county) , data=mn) 
sigma.y.hat <- sigma.hat(fit)$sigma$data
coef.hat <- as.matrix(coef(fit)$county["LAC QUI PARLE       ",])

x_new <- 1  #first floor
y_new <- rnorm(1000, coef.hat %*% c(1, x_new), sigma.y.hat  )
qplot(y_new)
```

What we got is what the log-radon measurements could possibly look like in Lac Qui Parle county. What about the radon measurements?

```{r message=F}
qplot(exp(y_new))
```

Let's get a 95% predictive interval

```{r}
quantile(y_new, c(0.025,  0.975))
quantile(exp(y_new), c(0.025,  0.975))
```

In our simulation, 95% of the data fell within the above interval.



Now, let's predict a new observation in a new county that we haven't seen.

```{r message=FALSE}
sigma_a <- 0.3487
sigma_b <- 0.3436
rho <- -0.34
Sigma <-  matrix(c(sigma_a^2           , rho*sigma_a*sigma_b, 
                   rho*sigma_a*sigma_b , sigma_b^2),
   nrow=2,       
   ncol=2,       
   byrow = TRUE)   

mu_a <- 1.4628
mu_b <- -0.6811
mu <- matrix(c(mu_a,
               mu_b),
               nrow=2,
               byrow=TRUE)
             
sigma_y <- 0.7462

random_params <- mvrnorm(n = 1000, mu, Sigma)


x_new = 1  # first floor

y_new <- rnorm(1000, mean=random_params %*% c(1, x_new), sd=sigma_y)
qplot(y_new)
```

Let's now construct a predictive interval again

```{r}
quantile(y_new, c(0.025,  0.975))
```


As you can see, the predictive interval is wider for an unknown county (of course).