# Divorce rates and their relationship with Marriage rate and Median Age Marriage


```R
# load data and copy
library(rethinking)
options(mc.cores = parallel::detectCores())
data(WaffleDivorce)
d <- WaffleDivorce
# standardize variables
d$A <- scale( d$MedianAgeMarriage )
d$D <- scale( d$Divorce )
```

    Loading required package: rstan
    
    Loading required package: StanHeaders
    
    Loading required package: ggplot2
    
    rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
    
    For execution on a local, multicore CPU with excess RAM we recommend calling
    options(mc.cores = parallel::detectCores()).
    To avoid recompilation of unchanged Stan programs, we recommend calling
    rstan_options(auto_write = TRUE)
    
    Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
    
    Loading required package: parallel
    
    rethinking (Version 2.13)
    
    
    Attaching package: 'rethinking'
    
    
    The following object is masked from 'package:stats':
    
        rstudent
    
    



```R
sd( d$MedianAgeMarriage )
```


1.24363030138808



```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 )
```

## Prior Predictive Simulation


```R
set.seed(10)
prior <- extract.prior( m5.1 )
mu <- link( m5.1 , post=prior , data=list( A=c(-2,2) ) )
plot( NULL , xlim=c(-2,2) , ylim=c(-2,2) )
for ( i in 1:50 ) lines( c(-2,2) , mu[i,] , col=col.alpha("black",0.4) )
```

## Divorce vs Median Age Marriage


```R
# compute percentile interval of mean
A_seq <- seq( from=-3 , to=3.2 , length.out=30 )
mu <- link( m5.1 , data=list(A=A_seq) )
mu.mean <- apply( mu , 2, mean )
mu.PI <- apply( mu , 2 , PI )
# plot it all
plot( D ~ A , data=d , col=rangi2 )
lines( A_seq , mu.mean , lwd=2 )
shade( mu.PI , A_seq )
```


```R
d$M <- scale( d$Marriage )
m5.2 <- quap(
alist(
D ~ dnorm( mu , sigma ) ,
mu <- a + bM * M ,
a ~ dnorm( 0 , 0.2 ) ,
bM ~ dnorm( 0 , 0.5 ) ,
sigma ~ dexp( 1 )
) , data = d )
```


```R
# compute percentile interval of mean
A_seq <- seq( from=-3 , to=3.2 , length.out=30 )
mu <- link( m5.2 , data=list(M=A_seq) )
mu.mean <- apply( mu , 2, mean )
mu.PI <- apply( mu , 2 , PI )
# plot it all
plot( D ~ M , data=d , col=rangi2 )
lines( A_seq , mu.mean , lwd=2 )
shade( mu.PI , A_seq )
```

# Directed Acyclic Graph


```R
library(dagitty)
dag5.1 <- dagitty( "dag {
A -> D
A -> M
M -> D
}")
coordinates(dag5.1) <- list( x=c(A=0,D=1,M=2) , y=c(A=0,D=1,M=0) )
drawdag( dag5.1 )
```

## Another DAG with a different causal implication


```R
DMA_dag2 <- dagitty('dag{ D <- A -> M }')
coordinates(DMA_dag2) <- list( x=c(A=0,D=1,M=2) , y=c(A=0,D=1,M=0) )
drawdag(DMA_dag2)
```


```R
impliedConditionalIndependencies( DMA_dag2 )
```


    D _||_ M | A


### The above relation translates to D is independent of A conditional on A

## The following is a multiple regression :-


```R
m5.3 <- quap(
alist(
D ~ dnorm( mu , sigma ) ,
mu <- a + bM*M + bA*A ,
a ~ dnorm( 0 , 0.2 ) ,
bM ~ dnorm( 0 , 0.5 ) ,
bA ~ dnorm( 0 , 0.5 ) ,
sigma ~ dexp( 1 )
) , data = d )
precis( m5.3 )

```


<table>
<caption>A precis: 4 × 4</caption>
<thead>
	<tr><th></th><th scope=col>mean</th><th scope=col>sd</th><th scope=col>5.5%</th><th scope=col>94.5%</th></tr>
	<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>
</thead>
<tbody>
	<tr><th scope=row>a</th><td> 6.145791e-06</td><td>0.09707503</td><td>-0.1551385</td><td> 0.1551508</td></tr>
	<tr><th scope=row>bM</th><td>-6.539899e-02</td><td>0.15077111</td><td>-0.3063603</td><td> 0.1755624</td></tr>
	<tr><th scope=row>bA</th><td>-6.135264e-01</td><td>0.15098166</td><td>-0.8548242</td><td>-0.3722285</td></tr>
	<tr><th scope=row>sigma</th><td> 7.851074e-01</td><td>0.07784077</td><td> 0.6607028</td><td> 0.9095120</td></tr>
</tbody>
</table>




```R
plot( coeftab(m5.1,m5.2,m5.3), par=c("bA","bM") )
```


```R
# call link without specifying new data
# so it uses original data
mu <- link( m5.3 )

```


```R
# summarize samples across cases
mu_mean <- apply( mu , 2 , mean )
mu_PI <- apply( mu , 2 , PI )
```


```R
mu_mean
```


<style>
.list-inline {list-style: none; margin:0; padding: 0}
.list-inline>li {display: inline-block}
.list-inline>li:not(:last-child)::after {content: "\00b7"; padding: 0 .5ex}
</style>
<ol class=list-inline><li>0.368354486649211</li><li>0.317846556171417</li><li>0.120172279734276</li><li>0.754594505305712</li><li>-0.35210242265142</li><li>0.114402384795662</li><li>-0.712026970344305</li><li>-0.322345564338635</li><li>-1.75748948800128</li><li>-0.118798644892876</li><li>0.0399074909436967</li><li>-0.501194958425354</li><li>1.30713399648199</li><li>-0.430038703852553</li><li>0.178068050755266</li><li>0.296693166231551</li><li>0.483538214776328</li><li>0.580402126397884</li><li>0.0657178960624549</li><li>-0.0585743662824406</li><li>-0.5847983864951</li><li>-1.13328867640735</li><li>-0.110195176519957</li><li>-0.0402545497768813</li><li>0.137379216480115</li><li>0.248008677498342</li><li>0.200437068524857</li><li>0.329386346048645</li><li>-0.310805774461407</li><li>-0.721743318476945</li><li>0.118451586059692</li><li>-1.10120331050511</li><li>0.167743888707763</li><li>0.256509397801258</li><li>-0.0677856485702234</li><li>0.750040238196824</li><li>0.04567738588231</li><li>-0.438034358310611</li><li>-0.971646219066463</li><li>-0.137726275313299</li><li>0.222198272379584</li><li>0.431412338693953</li><li>0.395277771527691</li><li>1.19245533419973</li><li>-0.354935996085724</li><li>-0.179022923503312</li><li>0.0519523466657839</li><li>0.481817521101744</li><li>-0.072947729593975</li><li>0.729896979946677</li></ol>




```R
mu_PI
```


<table>
<caption>A matrix: 2 × 50 of type dbl</caption>
<tbody>
	<tr><th scope=row>5%</th><td>0.1562013</td><td>-0.01750389</td><td>-0.04314372</td><td>0.4147012</td><td>-0.5525603</td><td>-0.1328497</td><td>-0.9793711</td><td>-0.648136479</td><td>-2.433501</td><td>-0.3490095</td><td>...</td><td>0.0448409</td><td>0.1612790</td><td>0.1981170</td><td>0.7121631</td><td>-0.5910848</td><td>-0.365989353</td><td>-0.1278724</td><td>0.2605910</td><td>-0.3022358</td><td>0.188137</td></tr>
	<tr><th scope=row>94%</th><td>0.5671654</td><td> 0.64324070</td><td> 0.28030955</td><td>1.0857921</td><td>-0.1537885</td><td> 0.3458035</td><td>-0.4343818</td><td>-0.006013426</td><td>-1.127705</td><td> 0.0894809</td><td>...</td><td>0.3966038</td><td>0.6773873</td><td>0.5886995</td><td>1.6518693</td><td>-0.1422689</td><td>-0.004714744</td><td> 0.2183203</td><td>0.6917975</td><td> 0.1341940</td><td>1.254396</td></tr>
</tbody>
</table>




```R
D_sim <- sim( m5.3 , n=1e4 )

```


```R
plot( mu_mean ~ d$D , col=rangi2 , ylim=range(mu_PI) ,
xlab="Observed divorce" , ylab="Predicted divorce" )
abline( a=0 , b=1 , lty=2 )
for ( i in 1:nrow(d) ) lines( rep(d$D[i],2) , mu_PI[,i] , col=rangi2 )
identify( x=d$D , y=mu_mean , labels=d$Loc )
```


```R
N <- 100 # number of cases
x_real <- rnorm( N ) # x_real as Gaussian with mean 0 and stddev 1
x_spur <- rnorm( N , x_real ) # x_spur as Gaussian with mean=x_real
y <- rnorm( N , x_real ) # y as Gaussian with mean=x_real
d <- data.frame(y,x_real,x_spur) # bind all together in data frame
```


```R
x_spur
```


<style>
.list-inline {list-style: none; margin:0; padding: 0}
.list-inline>li {display: inline-block}
.list-inline>li:not(:last-child)::after {content: "\00b7"; padding: 0 .5ex}
</style>
<ol class=list-inline><li>-1.35205556179838</li><li>0.692630018489859</li><li>2.38970488462669</li><li>-0.277785286445208</li><li>2.20878196395475</li><li>-0.623157211371598</li><li>-0.691301777675285</li><li>-0.00841657520744366</li><li>-2.1068966619116</li><li>-2.475596402338</li><li>3.77752924711212</li><li>0.0869765810510864</li><li>-2.63222491112067</li><li>-2.72039080437196</li><li>-0.663492280904649</li><li>1.92060795532126</li><li>0.335483979433785</li><li>-0.282782624992253</li><li>1.01393662342351</li><li>2.97442700466915</li><li>0.901401246933798</li><li>-2.05511179248847</li><li>0.765853777414379</li><li>-1.4979684583892</li><li>0.0813784389131222</li><li>-1.31675350397186</li><li>0.70254161891636</li><li>0.574591348151105</li><li>2.8984309960213</li><li>-1.55120387921632</li><li>0.624935822546091</li><li>0.588252047530261</li><li>1.32475039281867</li><li>-1.1329325267783</li><li>0.157429074311472</li><li>0.0186486883598361</li><li>-0.417081802417184</li><li>-1.72451659683524</li><li>1.93998990465857</li><li>0.243119622903495</li><li>-0.343638634450165</li><li>2.39354237507181</li><li>0.5500243128989</li><li>-1.18004489026227</li><li>-0.871716808627928</li><li>0.859613671896576</li><li>0.033162777666438</li><li>2.40866859364035</li><li>1.23765568960264</li><li>-1.31312755104094</li><li>0.0405348814554642</li><li>0.433098259672419</li><li>0.575605324917927</li><li>0.0515943799359837</li><li>-1.39389954009737</li><li>-0.239471665560266</li><li>-1.09359506382444</li><li>1.12313077723405</li><li>1.19603301999564</li><li>3.18237860849344</li><li>3.01756658487018</li><li>-0.101827794104095</li><li>0.00973837101934016</li><li>-0.89691510430296</li><li>0.799666659022323</li><li>-0.966017873812466</li><li>0.72400599076021</li><li>2.34177963098057</li><li>-0.168071738321592</li><li>-2.0217800588304</li><li>-1.93221396839229</li><li>3.3250880218069</li><li>-0.728398228232703</li><li>0.267828899677299</li><li>-0.656170322463872</li><li>-0.0455085821980248</li><li>-0.719921621345485</li><li>-1.16773318714732</li><li>1.77104102706098</li><li>-0.490194900479153</li><li>-0.219563777817322</li><li>-0.118072362754051</li><li>-0.503507180341146</li><li>1.27167119844878</li><li>-0.57045068299662</li><li>2.08256023802364</li><li>-1.66558305864451</li><li>1.6780100832135</li><li>0.285611914944682</li><li>-1.38461374984414</li><li>1.18848395567335</li><li>1.03807732050792</li><li>1.28997891579462</li><li>0.727116348329177</li><li>1.00115597456983</li><li>4.82655855587325</li><li>-0.412610227194</li><li>-1.96255662668122</li><li>1.58977310437138</li><li>0.869737378417309</li></ol>




```R
x_real
```


<style>
.list-inline {list-style: none; margin:0; padding: 0}
.list-inline>li {display: inline-block}
.list-inline>li:not(:last-child)::after {content: "\00b7"; padding: 0 .5ex}
</style>
<ol class=list-inline><li>0.46504871558418</li><li>0.104996326852726</li><li>0.591889277022188</li><li>0.456121018594855</li><li>1.13802392123417</li><li>-1.21721042748313</li><li>0.600311221961509</li><li>-1.10136482956525</li><li>-0.883278574282695</li><li>-0.969457867094617</li><li>2.21120726655929</li><li>0.781330749204854</li><li>-1.1628090726602</li><li>0.0164813378689372</li><li>-0.772072563770278</li><li>1.31261616585865</li><li>-1.30097706125507</li><li>1.89685450984775</li><li>1.55009918834605</li><li>1.51939375074762</li><li>-0.331432676279169</li><li>-1.66110778416189</li><li>0.712376193541223</li><li>-2.08962097413881</li><li>-0.052980535634609</li><li>-1.63824858304232</li><li>-0.352844033152804</li><li>0.690809926272784</li><li>1.3520587482415</li><li>-0.569337959990894</li><li>2.07566147817585</li><li>0.0401866561368483</li><li>0.0521689815121645</li><li>-0.957344116337265</li><li>0.000424660039592618</li><li>-0.574711312557669</li><li>-0.391885585632191</li><li>0.286070510211621</li><li>1.25198155374721</li><li>-0.429467160199921</li><li>-0.757008186470816</li><li>0.671992175500776</li><li>0.607393564444157</li><li>-0.0771225890925739</li><li>0.150965323857022</li><li>-0.452397772119638</li><li>-0.0634736510063597</li><li>-0.513526218325328</li><li>0.871929643061742</li><li>-0.963613924115809</li><li>0.354298983855878</li><li>-0.184439846755855</li><li>-0.143818257682809</li><li>1.38981428343578</li><li>-0.314723907597331</li><li>-1.36336930310008</li><li>0.558470264126142</li><li>-0.970500150066699</li><li>-0.00744609215997565</li><li>0.626770295873479</li><li>2.33889939043286</li><li>-0.176170200171509</li><li>-0.753746568899285</li><li>0.0674839915000366</li><li>-0.119438611501735</li><li>-0.742946582232449</li><li>0.306037345940029</li><li>1.45547127444296</li><li>0.4504564311188</li><li>-0.414867966712541</li><li>-1.2709710013827</li><li>1.49866852600129</li><li>0.351263269924126</li><li>-0.0141930404161029</li><li>-0.622948360782438</li><li>0.687860010292194</li><li>-0.618357566694332</li><li>-1.45181357454632</li><li>1.59663329619412</li><li>0.296462825072267</li><li>-0.407069732214744</li><li>0.617317395587624</li><li>-0.440971688864465</li><li>0.309480565407576</li><li>-0.404878083834447</li><li>-0.371848684746369</li><li>-1.32555213028749</li><li>1.52328907325857</li><li>0.355501079526269</li><li>0.568934108221817</li><li>0.680067387674436</li><li>0.219846401737705</li><li>1.34839467338858</li><li>0.0546604981055179</li><li>0.879076083934329</li><li>3.43648934510459</li><li>-0.726564047459416</li><li>-0.660313922161985</li><li>0.244668789049531</li><li>0.333072383646635</li></ol>




```R
d
```


<table>
<caption>A data.frame: 100 × 3</caption>
<thead>
	<tr><th scope=col>y</th><th scope=col>x_real</th><th scope=col>x_spur</th></tr>
	<tr><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>
</thead>
<tbody>
	<tr><td> 1.23837003</td><td> 0.46504872</td><td>-1.352055562</td></tr>
	<tr><td> 0.21014531</td><td> 0.10499633</td><td> 0.692630018</td></tr>
	<tr><td> 0.75867059</td><td> 0.59188928</td><td> 2.389704885</td></tr>
	<tr><td> 1.20829338</td><td> 0.45612102</td><td>-0.277785286</td></tr>
	<tr><td> 2.85487807</td><td> 1.13802392</td><td> 2.208781964</td></tr>
	<tr><td>-2.23129968</td><td>-1.21721043</td><td>-0.623157211</td></tr>
	<tr><td> 0.38395229</td><td> 0.60031122</td><td>-0.691301778</td></tr>
	<tr><td>-1.37574356</td><td>-1.10136483</td><td>-0.008416575</td></tr>
	<tr><td>-0.92992473</td><td>-0.88327857</td><td>-2.106896662</td></tr>
	<tr><td>-0.09839109</td><td>-0.96945787</td><td>-2.475596402</td></tr>
	<tr><td> 0.55906662</td><td> 2.21120727</td><td> 3.777529247</td></tr>
	<tr><td> 2.30905325</td><td> 0.78133075</td><td> 0.086976581</td></tr>
	<tr><td>-0.36630367</td><td>-1.16280907</td><td>-2.632224911</td></tr>
	<tr><td>-1.32368322</td><td> 0.01648134</td><td>-2.720390804</td></tr>
	<tr><td>-2.17579557</td><td>-0.77207256</td><td>-0.663492281</td></tr>
	<tr><td> 0.31392467</td><td> 1.31261617</td><td> 1.920607955</td></tr>
	<tr><td>-2.40565252</td><td>-1.30097706</td><td> 0.335483979</td></tr>
	<tr><td> 2.74501440</td><td> 1.89685451</td><td>-0.282782625</td></tr>
	<tr><td>-0.48731314</td><td> 1.55009919</td><td> 1.013936623</td></tr>
	<tr><td> 1.31195999</td><td> 1.51939375</td><td> 2.974427005</td></tr>
	<tr><td>-0.32518737</td><td>-0.33143268</td><td> 0.901401247</td></tr>
	<tr><td>-1.55476035</td><td>-1.66110778</td><td>-2.055111792</td></tr>
	<tr><td>-0.65982321</td><td> 0.71237619</td><td> 0.765853777</td></tr>
	<tr><td>-2.78895348</td><td>-2.08962097</td><td>-1.497968458</td></tr>
	<tr><td> 0.06102915</td><td>-0.05298054</td><td> 0.081378439</td></tr>
	<tr><td>-3.06322507</td><td>-1.63824858</td><td>-1.316753504</td></tr>
	<tr><td> 0.38972830</td><td>-0.35284403</td><td> 0.702541619</td></tr>
	<tr><td> 0.34871202</td><td> 0.69080993</td><td> 0.574591348</td></tr>
	<tr><td>-0.68836123</td><td> 1.35205875</td><td> 2.898430996</td></tr>
	<tr><td>-1.61809608</td><td>-0.56933796</td><td>-1.551203879</td></tr>
	<tr><td>...</td><td>...</td><td>...</td></tr>
	<tr><td> 0.06283928</td><td>-1.27097100</td><td>-1.93221397</td></tr>
	<tr><td> 2.24122861</td><td> 1.49866853</td><td> 3.32508802</td></tr>
	<tr><td>-0.90951450</td><td> 0.35126327</td><td>-0.72839823</td></tr>
	<tr><td>-1.18210173</td><td>-0.01419304</td><td> 0.26782890</td></tr>
	<tr><td>-2.73401509</td><td>-0.62294836</td><td>-0.65617032</td></tr>
	<tr><td> 1.54120743</td><td> 0.68786001</td><td>-0.04550858</td></tr>
	<tr><td> 0.55025982</td><td>-0.61835757</td><td>-0.71992162</td></tr>
	<tr><td>-1.62260770</td><td>-1.45181357</td><td>-1.16773319</td></tr>
	<tr><td> 2.02855006</td><td> 1.59663330</td><td> 1.77104103</td></tr>
	<tr><td> 0.99111728</td><td> 0.29646283</td><td>-0.49019490</td></tr>
	<tr><td>-1.74503457</td><td>-0.40706973</td><td>-0.21956378</td></tr>
	<tr><td> 0.91623296</td><td> 0.61731740</td><td>-0.11807236</td></tr>
	<tr><td> 0.88907842</td><td>-0.44097169</td><td>-0.50350718</td></tr>
	<tr><td> 2.17097253</td><td> 0.30948057</td><td> 1.27167120</td></tr>
	<tr><td> 0.55394043</td><td>-0.40487808</td><td>-0.57045068</td></tr>
	<tr><td>-1.15472825</td><td>-0.37184868</td><td> 2.08256024</td></tr>
	<tr><td>-0.16322237</td><td>-1.32555213</td><td>-1.66558306</td></tr>
	<tr><td> 1.61534860</td><td> 1.52328907</td><td> 1.67801008</td></tr>
	<tr><td> 0.32711912</td><td> 0.35550108</td><td> 0.28561191</td></tr>
	<tr><td>-0.38844307</td><td> 0.56893411</td><td>-1.38461375</td></tr>
	<tr><td> 1.24643241</td><td> 0.68006739</td><td> 1.18848396</td></tr>
	<tr><td>-0.26187694</td><td> 0.21984640</td><td> 1.03807732</td></tr>
	<tr><td> 1.72864891</td><td> 1.34839467</td><td> 1.28997892</td></tr>
	<tr><td> 0.72428338</td><td> 0.05466050</td><td> 0.72711635</td></tr>
	<tr><td>-0.03436937</td><td> 0.87907608</td><td> 1.00115597</td></tr>
	<tr><td> 3.40677569</td><td> 3.43648935</td><td> 4.82655856</td></tr>
	<tr><td>-0.45170983</td><td>-0.72656405</td><td>-0.41261023</td></tr>
	<tr><td> 1.68005195</td><td>-0.66031392</td><td>-1.96255663</td></tr>
	<tr><td> 0.37582248</td><td> 0.24466879</td><td> 1.58977310</td></tr>
	<tr><td>-0.16967498</td><td> 0.33307238</td><td> 0.86973738</td></tr>
</tbody>
</table>




```R
pairs(d)
```

# Counterfactual Plots

## Considering the below DAG we try to simulate the effect of A on D while also including the effect A has on M 


```R
drawdag( dag5.1 )
```


```R
data(WaffleDivorce)
d <- list()
d$A <- standardize( WaffleDivorce$MedianAgeMarriage )
d$D <- standardize( WaffleDivorce$Divorce )
d$M <- standardize( WaffleDivorce$Marriage )
m5.3_A <- quap(
alist(
## A -> D <- M
D ~ dnorm( mu , sigma ) ,
mu <- a + bM*M + bA*A ,
a ~ dnorm( 0 , 0.2 ) ,
bM ~ dnorm( 0 , 0.5 ) ,
bA ~ dnorm( 0 , 0.5 ) ,
sigma ~ dexp( 1 ),
## A -> M
M ~ dnorm( mu_M , sigma_M ),
mu_M <- aM + bAM*A,
aM ~ dnorm( 0 , 0.2 ),
bAM ~ dnorm( 0 , 0.5 ),
sigma_M ~ dexp( 1 )
) , data = d )
```


```R
A_seq <- seq( from=-2 , to=2 , length.out=30 )
```


```R
# prep data
sim_dat <- data.frame( A=A_seq )
```

## Given the data(A_seq) we simulate M first and later use this A_seq,M to generate D


```R
# simulate M and then D, using A_seq
s <- sim( m5.3_A , data=sim_dat , vars=c("M","D") )
```


```R
# display counterfactual predictions
plot( sim_dat$A , colMeans(s$D) , ylim=c(-2,2) , type="l" ,
xlab="manipulated A" , ylab="counterfactual D" )
shade( apply(s$D,2,PI) , sim_dat$A )
mtext( "Total counterfactual effect of A on D" )
```


```R
plot( sim_dat$A , colMeans(s$M) , ylim=c(-2,2) , type="l" ,
xlab="manipulated A" , ylab="counterfactual M" )
shade( apply(s$M,2,PI) , sim_dat$A )
mtext( "Total counterfactual effect of A on M" )
```

## The above counterfactual plots show that A has a direct influence on both A and M as assumed by our model

## What happends when we try to manipulate A without involving the causal change in M ?
1. Implies that we use A -> D and M -> D 
2. The code uses simulating **only D** assuming A = 0 (average State in US )


```R
# simulate D, using A_seq
sim_dat <- data.frame( M=seq(from=-2,to=2,length.out=30) , A=0 )
s <- sim( m5.3_A , data=sim_dat , vars="D")
```


```R
plot( sim_dat$M , colMeans(s) , ylim=c(-2,2) , type="l" ,
xlab="manipulated M" , ylab="counterfactual D" )
shade( apply(s,2,PI) , sim_dat$M )
mtext( "Total counterfactual effect of M on D" )
```

### The above plot shows that M doesnt have a significant affect on D

# Spurious Waffles 


```R
library(rethinking)
data(milk)
d <- milk
str(d)
```

    Loading required package: rstan
    
    Loading required package: StanHeaders
    
    Loading required package: ggplot2
    
    rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
    
    For execution on a local, multicore CPU with excess RAM we recommend calling
    options(mc.cores = parallel::detectCores()).
    To avoid recompilation of unchanged Stan programs, we recommend calling
    rstan_options(auto_write = TRUE)
    
    Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
    
    Loading required package: parallel
    
    rethinking (Version 2.13)
    
    
    Attaching package: 'rethinking'
    
    
    The following object is masked from 'package:stats':
    
        rstudent
    
    


    'data.frame':	29 obs. of  8 variables:
     $ clade         : Factor w/ 4 levels "Ape","New World Monkey",..: 4 4 4 4 4 2 2 2 2 2 ...
     $ species       : Factor w/ 29 levels "A palliata","Alouatta seniculus",..: 11 8 9 10 16 2 1 6 28 27 ...
     $ kcal.per.g    : num  0.49 0.51 0.46 0.48 0.6 0.47 0.56 0.89 0.91 0.92 ...
     $ perc.fat      : num  16.6 19.3 14.1 14.9 27.3 ...
     $ perc.protein  : num  15.4 16.9 16.9 13.2 19.5 ...
     $ perc.lactose  : num  68 63.8 69 71.9 53.2 ...
     $ mass          : num  1.95 2.09 2.51 1.62 2.19 5.25 5.37 2.51 0.71 0.68 ...
     $ neocortex.perc: num  55.2 NA NA NA NA ...



```R
d$K <- scale( d$kcal.per.g )
d$N <- scale( d$neocortex.perc )
d$M <- scale( log(d$mass) )
```


```R
m5.5_draft <- quap(
alist(
K ~ dnorm( mu , sigma ) ,
mu <- a + bN*N ,
a ~ dnorm( 0 , 1 ) ,
bN ~ dnorm( 0 , 1 ) ,
sigma ~ dexp( 1 )
) , data=d )
```


```R
d$neocortex.perc
```


<style>
.list-inline {list-style: none; margin:0; padding: 0}
.list-inline>li {display: inline-block}
.list-inline>li:not(:last-child)::after {content: "\00b7"; padding: 0 .5ex}
</style>
<ol class=list-inline><li>55.16</li><li>&lt;NA&gt;</li><li>&lt;NA&gt;</li><li>&lt;NA&gt;</li><li>&lt;NA&gt;</li><li>64.54</li><li>64.54</li><li>67.64</li><li>&lt;NA&gt;</li><li>68.85</li><li>58.85</li><li>61.69</li><li>60.32</li><li>&lt;NA&gt;</li><li>&lt;NA&gt;</li><li>69.97</li><li>&lt;NA&gt;</li><li>70.41</li><li>&lt;NA&gt;</li><li>73.4</li><li>&lt;NA&gt;</li><li>67.53</li><li>&lt;NA&gt;</li><li>71.26</li><li>72.6</li><li>&lt;NA&gt;</li><li>70.24</li><li>76.3</li><li>75.49</li></ol>




```R
dcc <- d[ complete.cases(d$K,d$N,d$M) , ]

```


```R
m5.5_draft <- quap(
alist(
K ~ dnorm( mu , sigma ) ,
mu <- a + bN*N ,
a ~ dnorm( 0 , 1 ) ,
bN ~ dnorm( 0 , 1 ) ,
sigma ~ dexp( 1 )
) , data=dcc )

```


```R
prior <- extract.prior( m5.5_draft )
xseq <- c(-2,2)
mu <- link( m5.5_draft , post=prior , data=list(N=xseq) )
plot( NULL , xlim=xseq , ylim=xseq )
for ( i in 1:50 ) lines( xseq , mu[i,] , col=col.alpha("black",0.3) )
```

## Impossible priors 

## Tightening the priors so that we have reasonable values


```R
m5.5 <- quap(
alist(
K ~ dnorm( mu , sigma ) ,
mu <- a + bN*N ,
a ~ dnorm( 0 , 0.2 ) ,
bN ~ dnorm( 0 , 0.5 ) ,
sigma ~ dexp( 1 )
) , data=dcc )
```


```R
prior <- extract.prior( m5.5 )
xseq <- c(-2,2)
mu <- link( m5.5_draft , post=prior , data=list(N=xseq) )
plot( NULL , xlim=xseq , ylim=xseq )
for ( i in 1:50 ) lines( xseq , mu[i,] , col=col.alpha("black",0.3) )
```


```R
precis( m5.5 )
```


<table>
<caption>A precis: 3 × 4</caption>
<thead>
	<tr><th></th><th scope=col>mean</th><th scope=col>sd</th><th scope=col>5.5%</th><th scope=col>94.5%</th></tr>
	<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>
</thead>
<tbody>
	<tr><th scope=row>a</th><td>0.03993997</td><td>0.1544908</td><td>-0.2069662</td><td>0.2868461</td></tr>
	<tr><th scope=row>bN</th><td>0.13323493</td><td>0.2237469</td><td>-0.2243559</td><td>0.4908258</td></tr>
	<tr><th scope=row>sigma</th><td>0.99982070</td><td>0.1647082</td><td> 0.7365852</td><td>1.2630562</td></tr>
</tbody>
</table>




```R
xseq <- seq( from=min(dcc$N)-0.15 , to=max(dcc$N)+0.15 , length.out=30 )
mu <- link( m5.5 , data=list(N=xseq) )
mu_mean <- apply(mu,2,mean)
mu_PI <- apply(mu,2,PI)
plot( K ~ N , data=dcc )
lines( xseq , mu_mean , lwd=2 )
shade( mu_PI , xseq )
```


```R
m5.6 <- quap(
alist(
K ~ dnorm( mu , sigma ) ,
mu <- a + bM*M ,
a ~ dnorm( 0 , 0.2 ) ,
bM ~ dnorm( 0 , 0.5 ) ,
sigma ~ dexp( 1 )
) , data=dcc )
precis(m5.6)
```


<table>
<caption>A precis: 3 × 4</caption>
<thead>
	<tr><th></th><th scope=col>mean</th><th scope=col>sd</th><th scope=col>5.5%</th><th scope=col>94.5%</th></tr>
	<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>
</thead>
<tbody>
	<tr><th scope=row>a</th><td> 0.04639191</td><td>0.1512784</td><td>-0.1953801</td><td>0.28816396</td></tr>
	<tr><th scope=row>bM</th><td>-0.28249661</td><td>0.1928751</td><td>-0.5907483</td><td>0.02575509</td></tr>
	<tr><th scope=row>sigma</th><td> 0.94923665</td><td>0.1570432</td><td> 0.6982513</td><td>1.20022199</td></tr>
</tbody>
</table>




```R
xseq <- seq( from=min(dcc$M)-0.15 , to=max(dcc$M)+0.15 , length.out=30 )
mu <- link( m5.6 , data=list(M=xseq) )
mu_mean <- apply(mu,2,mean)
mu_PI <- apply(mu,2,PI)
plot( K ~ M , data=dcc )
lines( xseq , mu_mean , lwd=2 )
shade( mu_PI , xseq )
```


```R
m5.7 <- quap(
alist(
K ~ dnorm( mu , sigma ) ,
mu <- a + bN*N + bM*M ,
a ~ dnorm( 0 , 0.2 ) ,
bN ~ dnorm( 0 , 0.5 ) ,
bM ~ dnorm( 0 , 0.5 ) ,
sigma ~ dexp( 1 )
) , data=dcc )
precis(m5.7)

```


<table>
<caption>A precis: 4 × 4</caption>
<thead>
	<tr><th></th><th scope=col>mean</th><th scope=col>sd</th><th scope=col>5.5%</th><th scope=col>94.5%</th></tr>
	<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>
</thead>
<tbody>
	<tr><th scope=row>a</th><td> 0.06799076</td><td>0.1340005</td><td>-0.1461680</td><td> 0.2821495</td></tr>
	<tr><th scope=row>bN</th><td> 0.67510201</td><td>0.2483044</td><td> 0.2782636</td><td> 1.0719405</td></tr>
	<tr><th scope=row>bM</th><td>-0.70298796</td><td>0.2207911</td><td>-1.0558547</td><td>-0.3501212</td></tr>
	<tr><th scope=row>sigma</th><td> 0.73803237</td><td>0.1324698</td><td> 0.5263200</td><td> 0.9497447</td></tr>
</tbody>
</table>




```R
plot( coeftab( m5.5 , m5.6 , m5.7 ) , pars=c("bM","bN") )
```

## Using both M and N have increased their influence on the outcome in a more farther direction


```R
 pairs( ~K + M + N ,
dcc )
```

## Counterfactual with N = 0


```R
xseq <- seq( from=min(dcc$M)-0.15 , to=max(dcc$M)+0.15 , length.out=30 )
mu <- link( m5.7 , data=data.frame( M=xseq , N=0 ) )
mu_mean <- apply(mu,2,mean)
mu_PI <- apply(mu,2,PI)
plot( NULL , xlim=range(dcc$M) , ylim=range(dcc$K) )
lines( xseq , mu_mean , lwd=2 )
shade( mu_PI , xseq )
```

## Counterfactual with M = 0


```R
xseq <- seq( from=min(dcc$N)-0.15 , to=max(dcc$N)+0.15 , length.out=30 )
mu <- link( m5.7 , data=data.frame( N=xseq , M=0 ) )
mu_mean <- apply(mu,2,mean)
mu_PI <- apply(mu,2,PI)
plot( NULL , xlim=range(dcc$N) , ylim=range(dcc$K) )
lines( xseq , mu_mean , lwd=2 )
shade( mu_PI , xseq )
```

## So it is clear that M and N when applied together lead to stronger effect 

# Index Variables

### Using dummy variables for categorical data doesnot help with the prior so we try to use Index variables


```R
data(Howell1)
d <- Howell1
str(d)
```

    'data.frame':	544 obs. of  4 variables:
     $ height: num  152 140 137 157 145 ...
     $ weight: num  47.8 36.5 31.9 53 41.3 ...
     $ age   : num  63 63 65 41 51 35 32 27 19 54 ...
     $ male  : int  1 0 0 1 0 1 0 1 0 1 ...



```R
mu_female <- rnorm(1e4,178,20)
mu_male <- rnorm(1e4,178,20) + rnorm(1e4,0,10)
precis( data.frame( mu_female , mu_male ) )
```


<table>
<caption>A precis: 2 × 5</caption>
<thead>
	<tr><th></th><th scope=col>mean</th><th scope=col>sd</th><th scope=col>5.5%</th><th scope=col>94.5%</th><th scope=col>histogram</th></tr>
	<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th></tr>
</thead>
<tbody>
	<tr><th scope=row>mu_female</th><td>177.6473</td><td>19.95584</td><td>145.6448</td><td>209.5282</td><td><span style=white-space:pre-wrap>&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2583&gt;&lt;U+2587&gt;&lt;U+2587&gt;&lt;U+2582&gt;&lt;U+2581&gt;&lt;U+2581&gt;                </span></td></tr>
	<tr><th scope=row>mu_male</th><td>178.2149</td><td>22.19615</td><td>142.7647</td><td>213.3100</td><td>&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2583&gt;&lt;U+2587&gt;&lt;U+2587&gt;&lt;U+2583&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;</td></tr>
</tbody>
</table>



### Prior predictive simulation clearly shows that male has more deviation in its values even though there is no prior reasoning to that

\begin{align}
h_i = Normal(\mu_i,\sigma)\\
\mu_i = \alpha_i + \beta_m*m_i 
\end{align}

since $\alpha_i$  represents the avg female height when $m_i$ = 0 and if it is 1 there are 2 parameters that determine male height causing more uncertainity even when there is no prior evidence


```R
d$sex <- ifelse( d$male==1 , 2 , 1 )
str( d$sex )
```

     num [1:544] 2 1 1 2 1 2 1 2 1 2 ...



```R
m5.8 <- quap(
alist(
height ~ dnorm( mu , sigma ) ,
mu <- a[sex] ,
a[sex] ~ dnorm( 178 , 20 ) ,
sigma ~ dunif( 0 , 50 )
) , data=d )

precis( m5.8 , depth=2 )
```


<table>
<caption>A precis: 3 × 4</caption>
<thead>
	<tr><th></th><th scope=col>mean</th><th scope=col>sd</th><th scope=col>5.5%</th><th scope=col>94.5%</th></tr>
	<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>
</thead>
<tbody>
	<tr><th scope=row>a[1]</th><td>134.91059</td><td>1.6069414</td><td>132.3424</td><td>137.47879</td></tr>
	<tr><th scope=row>a[2]</th><td>142.57873</td><td>1.6974815</td><td>139.8658</td><td>145.29164</td></tr>
	<tr><th scope=row>sigma</th><td> 27.31009</td><td>0.8280521</td><td> 25.9867</td><td> 28.63348</td></tr>
</tbody>
</table>




```R
post <- extract.samples(m5.8)
post$diff_fm <- post$a[,1] - post$a[,2]
precis( post , depth=2 )

```


<table>
<caption>A precis: 4 × 5</caption>
<thead>
	<tr><th></th><th scope=col>mean</th><th scope=col>sd</th><th scope=col>5.5%</th><th scope=col>94.5%</th><th scope=col>histogram</th></tr>
	<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th></tr>
</thead>
<tbody>
	<tr><th scope=row>sigma</th><td> 27.321169</td><td>0.819567</td><td> 26.01392</td><td> 28.619513</td><td>&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2583&gt;&lt;U+2587&gt;&lt;U+2587&gt;&lt;U+2587&gt;&lt;U+2583&gt;&lt;U+2582&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;</td></tr>
	<tr><th scope=row>a[1]</th><td>134.913639</td><td>1.619593</td><td>132.27227</td><td>137.503628</td><td><span style=white-space:pre-wrap>&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2582&gt;&lt;U+2585&gt;&lt;U+2587&gt;&lt;U+2587&gt;&lt;U+2585&gt;&lt;U+2582&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;        </span></td></tr>
	<tr><th scope=row>a[2]</th><td>142.546916</td><td>1.706710</td><td>139.82388</td><td>145.276203</td><td>&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2582&gt;&lt;U+2583&gt;&lt;U+2587&gt;&lt;U+2587&gt;&lt;U+2587&gt;&lt;U+2583&gt;&lt;U+2582&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;</td></tr>
	<tr><th scope=row>diff_fm</th><td> -7.633278</td><td>2.349692</td><td>-11.40781</td><td> -3.917746</td><td><span style=white-space:pre-wrap>&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2582&gt;&lt;U+2587&gt;&lt;U+2587&gt;&lt;U+2583&gt;&lt;U+2581&gt;&lt;U+2581&gt;&lt;U+2581&gt;                                </span></td></tr>
</tbody>
</table>



### So this variable  diff_fm  is called the Contrast


```R
data(milk)
d <- milk
unique(d$clade)
```


<style>
.list-inline {list-style: none; margin:0; padding: 0}
.list-inline>li {display: inline-block}
.list-inline>li:not(:last-child)::after {content: "\00b7"; padding: 0 .5ex}
</style>
<ol class=list-inline><li>Strepsirrhine</li><li>New World Monkey</li><li>Old World Monkey</li><li>Ape</li></ol>

<details>
	<summary style=display:list-item;cursor:pointer>
		<strong>Levels</strong>:
	</summary>
	<style>
	.list-inline {list-style: none; margin:0; padding: 0}
	.list-inline>li {display: inline-block}
	.list-inline>li:not(:last-child)::after {content: "\00b7"; padding: 0 .5ex}
	</style>
	<ol class=list-inline><li>'Ape'</li><li>'New World Monkey'</li><li>'Old World Monkey'</li><li>'Strepsirrhine'</li></ol>
</details>



```R
d$clade_id <- as.integer( d$clade )
d$clade_id
```


<style>
.list-inline {list-style: none; margin:0; padding: 0}
.list-inline>li {display: inline-block}
.list-inline>li:not(:last-child)::after {content: "\00b7"; padding: 0 .5ex}
</style>
<ol class=list-inline><li>4</li><li>4</li><li>4</li><li>4</li><li>4</li><li>2</li><li>2</li><li>2</li><li>2</li><li>2</li><li>2</li><li>2</li><li>2</li><li>2</li><li>3</li><li>3</li><li>3</li><li>3</li><li>3</li><li>3</li><li>1</li><li>1</li><li>1</li><li>1</li><li>1</li><li>1</li><li>1</li><li>1</li><li>1</li></ol>




```R
d$K <- scale( d$kcal.per.g )
m5.9 <- quap(
alist(
K ~ dnorm( mu , sigma ),
mu <- a[clade_id],
a[clade_id] ~ dnorm( 0 , 0.5 ),
sigma ~ dexp( 1 )
) , data=d )
labels <- paste( "a[" , 1:4 , "]:" , levels(d$clade) , sep="" )
plot( precis( m5.9 , depth=2 , pars="a" ) , labels=labels ,
xlab="expected kcal (std)" )

```

### The above example shows how this idea of Index Variables scales up 
And when you want to find the difference u need to simulate from the posterior distribution


```R

```

05af6c8344b00922449750019985506438594d3f
ipynb
Jupyter Notebook
The Many Variables & The Spurious Waffles.ipynb
GodEater8042/Statistical-Rethinking-Jupyter-R
ba305082b8fb24cefc43d02208de361e5adade3e
MIT
The Many Variables & The Spurious Waffles.ipynb
GodEater8042/Statistical-Rethinking-Jupyter-R
ba305082b8fb24cefc43d02208de361e5adade3e
MIT
The Many Variables & The Spurious Waffles.ipynb
GodEater8042/Statistical-Rethinking-Jupyter-R
ba305082b8fb24cefc43d02208de361e5adade3e
MIT
Qwen/Qwen-72B
1. YES
2. YES