Statistical Thinking 2023_v2.qmd

---
title: "Statistical Rethinking 2023"
author: "Luis Correia"
format: html
editor: visual
---

## Setup R

```{r}
library(rethinking)
library(tidyverse)
library(dplyr)
library(ggplot2)
library(rstan)
```

## Lecture 1

The content of this lecture should be reviewed and posted.

{{< pagebreak >}}

## Lecture 2

\<To be reviewed with lecture contents\>

A new model is born when we write:

$$
W\sim Binomial(N,p)
$$

with

$$
p\sim Unif(0,1)
$$

Simulation of Bayesian Experiment

```{r}
# definegrid
p_grid <-seq(from=0,to=1,length.out=20)

# defineprior
#prior <-rep(1,20)
prior <-ifelse(p_grid<0.5,0,1)

# computelikelihoodateachvalueingrid
likelihood <-dbinom(6,size=9,prob=p_grid)

# computeproductoflikelihoodandprior
unstd.posterior <-likelihood*prior

# standardizetheposterior,soitsumsto1
posterior <-unstd.posterior/sum(unstd.posterior)
```

```{r}
plot( p_grid,posterior,type="b",
xlab="probability ofwater",ylab="posteriorprobability")
mtext( "20points")
```

### Homework 1

Suppose the globe tossing data had turned into a 4-water and 11-land. Construct the posterior distribution.

```{r}
# 2.7 analyticalcalculation
W <-6
L <-3
curve( dbeta(x,W+1,L+1),from=0,to=1)

# quadraticapproximation
curve( dnorm(x,0.67,0.16),lty=2,add=TRUE)
```

Simulation of Bayesian Experiment

```{r}
# definegrid
p_grid <-seq(from=0,to=1,length.out=20)

# defineprior
#prior <-rep(1,20)
prior <-ifelse(p_grid<0.5,0,1)

# computelikelihoodateachvalueingrid
likelihood <-dbinom(6,size=9,prob=p_grid)

# computeproductoflikelihoodandprior
unstd.posterior <-likelihood*prior

# standardizetheposterior,soitsumsto1
posterior <-unstd.posterior/sum(unstd.posterior)
```

```{r}
plot( p_grid,posterior,type="b",
xlab="probability ofwater",ylab="posteriorprobability")
mtext( "20points")
```

### Homework 1

Suppose the globe tossing data had turned into a 4-water and 11-land. Construct the posterior distribution.

```{r}
# 2.7 analyticalcalculation
W <-6
L <-3
curve( dbeta(x,W+1,L+1),from=0,to=1)

# quadraticapproximation
curve( dnorm(x,0.67,0.16),lty=2,add=TRUE)
```

```{r}
p <- c(0, .25, .5, .75, 1)
model <- sapply(p, function (p, W, L) return (4*p)^W*(4-4*p)^L)
print(model)
```

```{r}
n_samples <-1000
p <-rep(NA,n_samples)
p[1] <-0.5
W <-6
L <-3
for (i in 2:n_samples){
  p_new <-rnorm(1,p[i-1],0.1)
  if (p_new < 0) p_new <- abs(p_new)
  if (p_new > 1) p_new <- 2-p_new
  q0 <-dbinom(W,W+L,p[i-1])
  q1 <-dbinom(W,W+L,p_new)
  p[i] <-ifelse(runif(1)<q1/q0,p_new,p[i-1])
}
```

```{r}
plot(density(p),xlim=c(0,1))
curve( dbeta(x, W+1, L+1 ),lty=2,add=TRUE)
```

```{r}
sample <- c('W','L','W','W','W','L','W','L','W')
W <- sum(sample=='W') # No. of W observed
L <- sum(sample=='L') # No. of L observed
p <- c(0,.25,.5,.75,1) # Proportions W
ways <- sapply(p, function (q) (q*4)^W*((1-q)*4)^L)
prob <- ways/sum(ways)
cbind(p, ways, prob)
```

```{r}
sim_globe <- function (p=.7, N = 9) {
  sample(c('W','L'), size = N, prob = c(p, 1-p), replace=TRUE)
}
sim_globe()
```

```{r}
replicate(sim_globe(p = .5, N = 9), n = 10)
```

```{r}
sum(sim_globe(p = .5, N = 1e4 )=='W')/1e4
```

```{r}
compute_posterior <- function (the_sample, poss = c(0,.25,.5,.75,1)) {
  W <- sum(the_sample=='W')  # No. of 'W'
  L <- sum(the_sample=='L')  # No. of 'L'
  ways <- sapply(p, function (q) (q*4)^W*((1-q)*4)^L)
  post <- ways/sum(ways)
  # bars <- sapply(post, function (q) make_bar(q))
  # data.frame(poss, ways, post=round(post,3), bars)
  df <- data.frame(poss, ways, post=round(post,3))
  ggplot(aes(x=poss, y=post), data=df)+
    geom_bar(stat="identity")
}
```

```{r}
the_sample <- sim_globe()

print(the_sample)

compute_posterior(the_sample)
```

Calculating the posteriors using the conjugate prior beta

```{r}

post_samples <- rbeta(1e3, 6+1, 3+1)
```

```{r}
dens(post_samples, lwd=4, col=2, xlab="proportion water", adj=.1)
curve(dbeta(x, 6+1, 3+1), add=TRUE, lty=2, lwd=3)
```

```{r}
w <-6;n<-9;
p_grid <-seq(from=0,to=1,length.out=100)
posterior <-dbinom(w,n,p_grid)*dunif(p_grid,0,1)
posterior <-posterior/sum(posterior)
plot(posterior,col=2, type="l")
```

## Lecture 3

### Workflow

> 1.  State a clear Question
> 2.  State your causal Assumptions
> 3.  Use the causal sketch to define the Generative Model
> 4.  Use the Generative Model to build an Estimator
> 5.  Profit

Using the Howell database, we will construct a model to analyse the causal relationship between Height and Weight.

A new model is born when we write:

$$
W = \beta H + U
$$

```{r}

sim_weight <- function (H, b, sd) {
  U <-  rnorm(length(H),0,sd)
  W <- b*H + U
  return(W)
}
```

```{r}
H <-  runif(200, min=130, max=170)
W <- sim_weight(H, b=.5, sd=5)
plot(W ~ H, col=2, lwd=3)
```

Our model can be written as follows

$$
W_i=\beta H_i+U_i
$$

$$
H_i\sim Uniform(130, 170)
$$

$$
U_i\sim Normal(0,\sigma^2)
$$

and our estimator can then be written as

$$
E(W_i|H_i)=\alpha+\beta H_i
$$

The posterior distribution is now written as\
$$
Pr(\alpha,\beta,\sigma^2|H_i,W_i)\propto Pr(W_i|H_i,\alpha,\beta,\sigma^2)Pr(\alpha,\beta,\sigma^2)
$$

Now our Bayesian Model is defined as follows:

$$
W_i\sim Normal(\mu_i,\sigma^2)
$$

$$
\mu_i = \alpha+\beta H_i
$$

```{r}

data(Howell1)
d <-Howell1[Howell1$age>=18,]
```

```{r}
#Similate a sample of 10 people
set.seed(93)
H <- runif(10,130,170)
W <- sim_weight(H,b=0.5, sd=5)

#Run Model
m3.1 <- quap(
  alist(
    W ~ dnorm(mu,sigma),
    mu <- a + b*H,
    a ~ dnorm(0,10),
    b ~ dunif(0,1),
    sigma~dunif(0,10)
),data=list(W=W, H=H))

#Summary
precis(m3.1)
```

```{r}
n <- 1e3
a <- rnorm(n,0,10)
b <- runif(n,0,1)
plot(NULL,xlim=c(130,170),ylim=c(50,90),
     xlab="height(cm)", ylab="weight(kg)",)
for (j in 1:50 ) abline(a=a[j], b=b[j], lwd=2, col=2)
```

```{r}
dat <- list(W=d$weight, H=d$height)
#k <- 10
#dat <- list(W=d$weight[1:k], H=d$height[1:k])

#Run Model
m3.2 <- quap(
  alist(
    W ~ dnorm(mu,sigma),
    mu <- a + b*H,
    a ~ dnorm(0,10),
    b ~ dunif(0,1),
    sigma~dunif(0,10)
),data=dat)
precis(m3.2)
post <- extract.samples(m3.2)
plot(dat$H,dat$W,col=2,lwd=3,
     xlim=c(min(dat$H),max(dat$H)), ylim=c(min(30,dat$w),max(dat$W)),
     xlab="height(cm)", ylab="weight(kg)")
for (j in 1:20)
  abline(a=post$a[j],b=post$b[j], lty=1)
```

```{r}
post <- extract.samples(m3.2)
pairs(post)
head(post)
plot(d$height,d$weight,col=2,lwd=3,
     xlab="height(cm)", ylab="weight(kg)")
for (j in 1:20)
  abline(a=post$a[j],b=post$b[j], lty=1)

height_seq <- seq(130,190,len=20)
W_postpred <- sim(m3.2,data=list(H=height_seq))
W_PI <- apply(W_postpred,2,PI)
lines(height_seq,W_PI[1,],lty=2,lwd=2)
lines(height_seq,W_PI[2,],lty=2,lwd=2)
```

From book (pp.51)

```{r}
pos_true <- .95 # likelihood
pos <- .01 # prior
true_pos <- (pos_true*pos)/(pos_true*pos+(1-pos_true)*(1-pos))
true_pos  # posterior
```

```{r}
p <- seq(from=.01, to=1, by=0.01); ptr <- .95
true_pos <- sapply(p, function (q) (ptr*q)/(ptr*q+(1-ptr)*(1-q)))

plot(true_pos, p, col=2, lwd=3)
```

### Intermezzo - Central Limit Theorem convergence

```{r}
n=20
logYs = -9.89
A = sum(logYs)-1
B = 2*n-1+sum(logYs)
C = n

theta_hat = (-B-sqrt(B^2-4*A*C))/(2*A) ##Quadratic formula

theta_vec = seq(0, 6, length=500)
post_prop = theta_vec^n*(theta_vec+1)^n*exp(-theta_vec)*(exp(sum(logYs)))^theta_vec

# pdf(file = "~/Dropbox/pubh8442/2015/Asymptotic_Posterior_Approximation_thetaDens.pdf", width = 7, height=4)

plot(theta_vec,post_prop,xlim = c(0,6),  ylab = "Unormalized density", xlab = "theta", type = 'l', lwd = 2, main = "")
abline(v=theta_hat, col = 'red')
#dev.off()

var = (1/20)*(theta_hat^2*(1+theta_hat)^2)/(theta_hat^2+(theta_hat+1)^2)
post_approx = dnorm(theta_vec, theta_hat,sqrt(var))

#pdf(file = "~/Dropbox/pubh8442/2015/Asymptotic_Posterior_Approximation_NormApp##.pdf", width = 7, height=4)

plot(theta_vec,post_approx,xlim = c(0,6),  ylab = "Density", xlab = "theta", type = 'l', lwd = 2, main = "")
abline(v=theta_hat, col = 'red')
#dev.off()
```

#### Difference between Percentile and C.I. for Mean

```{r}
n <- 100

nsurveys <- 200

meanSurveys <- vector()
sdSurveys <- vector()
SSampl <- vector()
percs <- vector()


for (i in 1:nsurveys) {
  S <- rnorm(n, mean=2, sd=5)
  mu_hat <- mean(S)
  sd_hat <- sd(S)
  SSampl <- c(SSampl, list(S))
  percs <- c(percs,list(quantile(S,probs=c(.025,.975))))
  meanSurveys <- c(meanSurveys, mu_hat)  # append mu_hat
  sdSurveys <- c(sdSurveys, sd_hat)   # append sd_hat
}

# Calculating MOE
MOE <- qnorm(.975)*sdSurveys/sqrt(n)

dataSampl <- data.frame (
  mu_hat = meanSurveys,
  sd_hat = sdSurveys,
  MOE = MOE,
  lower = meanSurveys-qnorm(.975)*sdSurveys/sqrt(n),
  upper = meanSurveys+qnorm(.975)*sdSurveys/sqrt(n)
)

```

```{r, fig.height = 3.5, fig.width = 6.5}
set.seed(42)
SSlice <- sample(1:nsurveys,100)
dataSampl %>%
  slice(SSlice) %>% # select 100 samples
  ggplot(aes(x = c(1:100))) + 
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.1) + 
  geom_point(aes(y=mu_hat)) + 
  xlab("Sample of Approval Estimate") + ylab("95% I.C.") + 
  theme_bw()
```

```{r}
i <- 5
L <- SSampl[SSlice[i]]
# quantile(L[[1]],probs=c(.025,.975))
percs[SSlice][i-1]
percs[SSlice][i]
```

## Lecture 4

### Workflow

```{r}
data(Howell1)
# d <-Howell1[Howell1$age>=18,]
d <-Howell1
```

```{r}
d |> 
  ggplot(aes(x=height, y=weight, color=factor(male)))+
  geom_point()+
  theme_bw()

d |> 
  ggplot(aes(x=weight))+
  geom_density(aes(color=factor(male)))+
  theme_bw()

d |> 
  ggplot(aes(x=height))+
  geom_density(aes(color=factor(male)))+
  theme_bw()
```

```{r}
# S = 1: female, S=2: male
sim_HW <- function(S, b, a) {
  N <- length(S)
  H <- ifelse(S==1,150,160)+rnorm(N,0,5)
  W <- a[S] + b[S]*H + rnorm(N, 0, 5)
  data.frame(S,H,W)
}
```

Simulating our data

```{r}
S <- rbern(100)+1
dat <- sim_HW(S, b=c(.5,.6), a=c(0,0))
head(dat)

dat |> 
  ggplot(aes(x=H, y=W, color=factor(S)))+
  geom_point()+
  theme_bw()
```

### Testing the Model

```{r}
# female sample
S <- rep(1,100)
simF <- sim_HW(S, b=c(.5,.6), a=c(0,0))

# male sample
S <- rep(2,100)
simM <- sim_HW(S, b=c(.5,.6), a=c(0,0))

# effect of sex (male-female)
mean(simM$W-simF$W)

```

### Determine the Estimator of our model

```{r}
S <- rbern(100)+1
dat <- sim_HW(S, b=c(.5,.6), a=c(0,0))

# estimate posterior
m_SW <- quap(
  alist(
    W~dnorm(mu,sigma),
    mu <- a[S],
    a[S] ~ dnorm(60,10),
    sigma ~ dunif(0,10)
  ), data=dat
)
precis(m_SW,depth=2)
```

### Running our model Over the Data

```{r}
# recoding
dat <- list(
  W = d$weight,
  S = d$male+1
)
# estimate posterior
m_SW <- quap(
  alist(
    W~dnorm(mu,sigma),
    mu <- a[S],
    a[S] ~ dnorm(60,10),
    sigma ~ dunif(0,10)
  ), data=dat
)
precis(m_SW,depth=2)
```

```{r}
post <- extract.samples(m_SW)
dens(post$a[,1],xlim=c(35, 55),lwd=3,col=2,xlab="posterior mean weight(kg")
dens(post$a[,2], lwd=3,col=4,add=TRUE)
```

Simulating the weights using the posteriors obtrained

```{r}
W1 <- rnorm(1000,post$a[,1], post$sigma)
W2 <- rnorm(1000,post$a[,2], post$sigma)
dens(W1,xlim=c(20,70),ylim=c(0,.085),lwd=3, col=2)
dens(W2,lwd=3, col=4, add=TRUE)

```

### Computing the Contrast

Now we need to compute the **contrast distribution**, the difference between the categories

```{r}
mu_contrast <- post$a[,2]-post$a[,1]
dens(mu_contrast,xlim=c(3,10),lwd=3, col=1,xlab="posterior mean weight contrast(kg)")
```

```{r}
W1 <- rnorm(1000,post$a[,1], post$sigma)
W2 <- rnorm(1000,post$a[,2], post$sigma)

# conntrast
W_contrast <- W2 - W1
dens(W_contrast,xlim=c(-25,35),lwd=3, col=1, 
     xlab="posterior weight contrast (kg)")

W_up <- sum(W_contrast>0)/1000
W_down <- sum(W_contrast<0)/1000
cat("W_up:",W_up," and W_down:",W_down)

df <- data.frame (
  W1 = W1,
  W2 = W2,
  W_contrast = W_contrast,
  W_flag = factor(W_contrast>W_up)
)

ggplot(df, aes(x=W_contrast,fill=W_flag))+
  geom_density()+
  theme_bw()
```

Now, what is the direct effect of sex over weight?

```{r, message=FALSE}
S <- rbern(100)+1
dat <- sim_HW(S, b=c(.5,.5), a=c(0,10))

# Indirect effect (b) is the same for men and women
dat |> 
  ggplot(aes(x=H, y=W, color=factor(S)))+
  geom_point()+
  geom_smooth(method = lm)+
  # facet_wrap(~S)+
  theme_bw()
```

```{r, message=FALSE}
dat <- list(
  W = d$weight,
  H = d$height,
  Hbar = mean(d$height),
  S = d$male+1 )

m_SHW <- quap(
  alist(
    W ~ dnorm(mu,sigma),
    mu <- a[S]+b[S]*(H-Hbar),
    a[S] ~ dnorm(60,10),
    b[S] ~ dunif(0,1),
    sigma ~ dunif(0,10)
  ), data=dat
)
data.frame(dat) |> 
  ggplot(aes(x=H, y=W, color=factor(S)))+
  geom_point()+
  geom_smooth(method = lm)+
  # facet_wrap(~S)+
  theme_bw()

precis(m_SHW,depth=2)
```

Note that regression line's slopes are almost the same for both male and female.

```{r}
post <- extract.samples(m_SHW)
dens(post$b[,1],xlim=c(.3, .9),ylim=c(0,8),
     lwd=3,col=2,xlab="posterior mean weight(kg)")
dens(post$b[,2],lwd=3,col=4,add=TRUE)
```

```{r}

```

### Polynomials

```{r}
plot(height~weight,d)
```

```{r}
d$weight_s <-(d$weight-mean(d$weight))/sd(d$weight)
d$weight_s2 <-d$weight_s^2
m4.5 <-quap(
  alist(
    height ~dnorm(mu,sigma),
    mu <-a+b1*weight_s+b2*weight_s2,
    a ~dnorm(178,20),
    b1 ~dlnorm(0,1),
    b2 ~dnorm(0,1),
    sigma ~dunif(0,50)
  ) ,data=d)
```

```{r}
precis(m4.5)
```

```{r}
# Not working - To be Reviewed
stan_m4.5 <- "
data {
  int<lower = 1> N;
  vector[N] height;
  vector[N] weight_s;
  vector[N] weight_s2;
}
parameters {
  real a;
  real b1;
  real b2;
  real <lower=0> sigma;
}

transformed parameters {
  vector[N] mu;
  
  for (i in 1:N)
    mu[i] = a + b1*weight_s[i] + b2*weight_s2[i];
}

model {
  // priors
  a ~ normal(178,20);
  b1 ~ normal(0,1);
  b2 ~ normal(0,1);
  sigma ~ uniform(0,50);
  
  // model  
  height ~normal(mu,sigma);
}"
data_m4.5  <- list(
  N = nrow(d),
  height = d$height,
  weight_s = d$weight_s,
  weight_s2 = d$weight_s2
)

#fit_m4.5 <- stan(data = data_m4.5, model_code = stan_m4.5, 
#                 chains = 3, iter = 1000)

```

```{r}
# Same Model in STAN
stan_m4.6 <- "
data {
  int<lower = 1> N;
  vector[N] height;
  vector[N] weight_s;
  vector[N] weight_s2;
}

transformed data{
  real bar_x1;
  real x1_sd;
  vector[N] x1_std;
  real y_sd;

  bar_x1 = mean(weight_s);
  x1_sd = sd(weight_s);
  x1_std = (weight_s - bar_x1)/x1_sd; // centered and scaled

  y_sd = sd(height);
}

parameters {
  real alpha_std;
  real beta1_std;
  real beta3_std;
  real <lower=0> sigma;
}

transformed parameters {
  vector[N] mu;
  
  for (i in 1:N)
    mu[i] = alpha_std + beta1_std*weight_s[i] + 
        beta3_std*weight_s2[i];
}

model {
  // priors
  alpha_std ~ normal(0,20);
  beta1_std ~ normal(0,2);
  beta3_std ~ normal(0,2);
  //sigma ~ normal(5,1.5);
  sigma ~ exponential(1/y_sd);  
  
  // model  
  height ~normal(mu,sigma);
}

generated quantities {
  real alpha;
  real beta1;
  real beta3;

  alpha = alpha_std - (beta1_std*bar_x1)/x1_sd 
      + (beta3_std*bar_x1^2)/x1_sd^2;

  beta1 = beta1_std/x1_sd - 2*(beta3_std*bar_x1)/x1_sd^2;

  beta3 = beta3_std/x1_sd^2;
}"

data_m4.6  <- list(
  N = nrow(d),
  height = d$height,
  weight_s = d$weight_s,
  weight_s2 = d$weight_s2
)

fit_m4.6 <- stan(data = data_m4.6, 
                 model_code = stan_m4.6,
                 chains = 4,
                 iter = 2000)

```

```{r}
print(fit_m4.6, pars = c("alpha","beta1", "beta3", "sigma"))
print(fit_m4.6, pars = c("alpha_std","beta1_std", "beta3_std", "sigma"))
precis(m4.5)
```

```{r}
traceplot(fit_m4.6, pars = "sigma")
```

```{r}
sigma_sample <- extract(fit_m4.6)[["sigma"]]
hist(sigma_sample)
```

### Cherry Blossoms example (pp.114) - Splines

```{r}
data(cherry_blossoms)
d <-cherry_blossoms
d2 <- d[ complete.cases(d$temp),]
precis(d2)
```

```{r}
?cherry_blossoms
d2 <-d[complete.cases(d$doy),]#completecasesondoy
num_knots <-15
knot_list <-quantile(d2$year,probs=seq(0,1,length.out=num_knots))
```

Building B (Block Variables)

```{r}
library(splines)
B <-bs(d2$year,
       knots=knot_list[-c(1,num_knots)] ,
       degree=3 ,intercept=TRUE)
```

```{r}
plot( NULL,xlim=range(d2$year),ylim=c(0,1),xlab="year",ylab="basis")
for (i in 1:ncol(B))lines(d2$year,B[,i])
```

```{r}
m4.7 <-quap(
  alist(
    D ~dnorm(mu,sigma),
    mu <-a+B%*%w,
    a ~dnorm(100,10),
    w ~dnorm(0,10),
    sigma ~dexp(1)
  ), data=list(D=d2$doy,B=B),
  start=list( w=rep(0,ncol(B))))
```

```{r}
precis(m4.7)
```

```{r}
post <-extract.samples(m4.7)
w <-apply(post$w,2,mean)
#plot( NULL,xlim=range(d2$year),ylim=c(-6,6),
#      xlab="year" ,ylab="basis*weight")
plot( d2$year,d2$doy-m4.7@optim$par[18],col=col.alpha(rangi2,0.3),pch=16,
      xlim=range(d2$year),ylim=c(-26,26),
      xlab="year" ,ylab="basis*weight")
for (i in 1:ncol(B))lines(d2$year,w[i]*B[,i])
```

```{r}
mu <-link(m4.7)
mu_PI <-apply(mu,2,PI,0.97)
plot( d2$year,d2$doy,col=col.alpha(rangi2,0.3),pch=16)
shade( mu_PI,d2$year,col=col.alpha("black",0.5))
```

```{r}
m4.7alt <-quap(
  alist(
    D ~dnorm(mu,sigma),
    mu <-a+sapply(1:827,function(i)sum(B[i,]*w)),
    a ~dnorm(100,1),
    w ~dnorm(0,10),
    sigma ~dexp(1)
  ),
  data=list( D=d2$doy,B=B),
  start=list( w=rep(0,ncol(B))))
```

Now using the B-Splines with the Weight x Height Data, we have the following:

```{r}
data(Howell1)
# d <-Howell1[Howell1$age>=18,]
d <-Howell1
```

## Lecture 5

```{r}
n <- 1000
Z <- rbern(n, .5)
X <- rbern(n, (1-Z)*.1+Z*.9)
Y <- rbern(n, (1-Z)*.1+Z*.9)
```

```{r}
cor(X,Y)
```

```{r}
cols <- c(4,2)

N <- 300
Z <- rbern(N)
X <- rnorm(N,2*Z-1)
Y <- rnorm(N,2*Z-1)

plot(X,Y, col=cols[Z+1], lwd=3)

abline(lm(Y[Z==1]~X[Z==1]), col=2, lwd=3)

abline(lm(Y[Z==0]~X[Z==0]), col=4, lwd=3)

abline(lm(Y~X), lwd=3)
```

Divorce Database and Waffle Houses

```{r}
data(WaffleDivorce)

d <-WaffleDivorce

# standardizevariables
d$D <-standardize(d$Divorce)
d$M <-standardize(d$Marriage)
d$A <-standardize(d$MedianAgeMarriage)

dat <- list(
  D = d$D,
  M = d$M,
  A = d$A
)
```

```{r}
m_DMA <- quap(
  alist(
    D ~ dnorm(mu, sigma),
    mu <- a + bM*M + bA*A,
    a ~dnorm(0,.2),
    bM~dnorm(0,.5),
    bA~dnorm(0,.5),
    sigma~dexp(1)
  ), data=dat )
```

```{r}
m5.1 <-quap(
  alist(
    D ~dnorm(mu,sigma),
    mu <-a+bA*A,
    a ~dnorm(0,0.2),
    bA ~dnorm(0,0.5),
    sigma ~dexp(1)
  ) ,data=d)
```

```{r}
precis(m_DMA)
plot(precis(m_DMA))
```

```{r}
precis(m5.1)
```