<jupyter_start><jupyter_text>Rethinking Statistics course in NumPyro - Week 7 Lecture 13: Monsters & Mixtures (Poisson GLMs, survival, zero-inflation)- [Video](https://www.youtube.com/watch?v=p7g-CgGCS34)- [Slides](https://speakerdeck.com/rmcelreath/l13-statistical-rethinking-winter-2019)Lecture 14: Ordered Categories, Left & Right- [Video](https://www.youtube.com/watch?v=zA3Jxv8LOrA)- [Slides](https://speakerdeck.com/rmcelreath/l14-statistical-rethinking-winter-2019)[Proposed problems](https://github.com/gbosquechacon/statrethinking_winter2019/blob/master/homework/week07.pdf) and [solutions in R](https://github.com/gbosquechacon/statrethinking_winter2019/blob/master/homework/week07_solutions.pdf) for the exercises of the week.<jupyter_code>import pandas as pd
import numpy as np

from jax import random
import jax.numpy as jnp
import numpyro
import numpyro.distributions as dist
from numpyro.distributions.transforms import OrderedTransform
from numpyro.infer import NUTS, MCMC, Predictive
%load_ext watermark
%watermark -n -u -v -iv -w
rng_key = random.PRNGKey(0)<jupyter_output><empty_output><jupyter_text>Exercise 1 > In the Trolley data, `Trolley`, we saw how education level (modeled as an ordered category) is associated with responses. Is this association causal? One plausible confound is that education is also associated with age, through a causal process: People are older when they finish school than when they begin it.> Reconsider the `Trolley` data in this light. Draw a DAG that represents hypothetical causal relationships among response, education, and age. Which statical model or models do you need to evaluate the causal influence of education on responses? Fit these models to the trolley data. What do you conclude about the causal relationships among these three variables? This is my DAG:Age could influence both response `R` and education `E`. It could influence `R`, because people at different ages could have different attitudes. Age could influence education, because the longer you have lived, the more education you could have completed (up to a point). It's like the age causing marriage example from earlier in the course.To evaluate the causal influence of `E` on `R`, we need to block the back-door from `E` through `A` to `R`. So we need a model that conditions on both `E` and `A`. Then the estimate for `E` should be the causal influence of `E`.<jupyter_code>d = pd.read_csv('./dat/Trolley.csv', header=0, sep=';')
elvl = d['edu'].unique()
idx = [7 , 0 , 6 , 4 , 2 , 1, 3, 5]
cat = pd.Categorical(d.edu, categories=list(elvl[idx]), ordered=True)
d['edu_cat'] = pd.factorize(cat, sort=True)[0].astype('int')
d['age_std'] = (d.age - d.age.mean())/d.age.std()
d.tail(3)
# https://fehiepsi.github.io/rethinking-numpyro/12-monsters-and-mixtures.html

def model(action, contact, intention, age_std, edu_cat, response=None):
    # priors
    cutpoints = numpyro.sample('cutpoints', 
                               dist.TransformedDistribution( 
                                   dist.Normal(0, 1.5).expand([6]), 
                                   OrderedTransform()
                                   ),
                               )
    norm = lambda label: numpyro.sample(label, dist.Normal(0,0.5))
    bA = norm('bA')
    bC = norm('bC')
    bI = norm('bI')
    bIA = norm('bIA')
    bIC = norm('bIC')
    bAge = norm('bAge')
    bEdu = norm('bEdu')
    edu_coefs = numpyro.sample('edu_coefs', dist.Dirichlet(jnp.repeat(2,7)))
    # likelihood
    delta_j = jnp.pad(edu_coefs, (1, 0))
    delta_E = jnp.sum(jnp.where(jnp.arange(8) <= edu_cat[..., None], delta_j, 0), -1)
    BI = bI + bIA*action + bIC*contact
    phi = bA*action + BI*intention + bC*contact + bAge*age_std + bEdu*delta_E
    response_hat = numpyro.sample('response_hat', dist.OrderedLogistic(phi, cutpoints), obs=response-1)
kernel = NUTS(model)
mcmc = MCMC(kernel, num_warmup=500, num_samples=500, num_chains=4, chain_method='sequential')
dat = {k:v.to_numpy() for k,v in d[['action', 'contact', 'intention', 'age_std', 'edu_cat', 'response']].items()}
mcmc.run(rng_key, **dat)
mcmc.print_summary()<jupyter_output>sample: 100%|██████████| 1000/1000 [01:14<00:00, 13.41it/s, 31 steps of size 1.41e-01. acc. prob=0.94]
sample: 100%|██████████| 1000/1000 [00:53<00:00, 18.57it/s, 31 steps of size 1.74e-01. acc. prob=0.92]
sample: 100%|██████████| 1000/1000 [00:48<00:00, 20.64it/s, 31 steps of size 1.42e-01. acc. prob=0.94]
sample: 100%|██████████| 1000/1000 [00:53<00:00, 18.66it/s, 63 steps of size 9.42e-02. acc. prob=0.95]<jupyter_text>You may recall from the chapter that education has a negative effect in the model without age. Now that we include age, education has a positive influence (with some overlap with zero). So age has indeed soaked up some of the previous influence assigned to education. The back-door may be real.I'd summarize this model, assuming this DAG is true, as saying that age causes people to give slightly lower responses. This could be a cohort effect, and not a causal influence of age. Either way, it is small. Education seems to cause higher responses (more approval). This suggests that education trains people to see some or all of the features A, I, C as more permissible. A model that interacted education with each might shed more light on things. Remember: A DAG doesn't say whether you need an interaction effect or not. That is a separate problem.  Exercise 2 > Consider one more variable in the Trolley data: Gender. Suppose that gender might influence education as well as response directly. Draw the DAG now that includes response, education, age, and gender.> Using only the DAG, is it possible that the inferences from Problem 1 are confounded by gender? If so, define any additional models you need to infer the causal influence of education on response. What do you conclude? This is my DAG:Here's the model we need, which includes education, age, and gender (female dummy variable):<jupyter_code>d = pd.read_csv('./dat/Trolley.csv', header=0, sep=';')
elvl = d['edu'].unique()
idx = [7 , 0 , 6 , 4 , 2 , 1, 3, 5]
cat = pd.Categorical(d.edu, categories=list(elvl[idx]), ordered=True)
d['edu_cat'] = pd.factorize(cat, sort=True)[0].astype('int')
d['age_std'] = (d.age - d.age.mean())/d.age.std()
d.tail(3)
def model(action, contact, intention, age_std, edu_cat, male, response=None):
    # priors
    cutpoints = numpyro.sample('cutpoints', 
                               dist.TransformedDistribution( 
                                   dist.Normal(0, 1.5).expand([6]), 
                                   OrderedTransform()
                                   ),
                               )
    norm = lambda label: numpyro.sample(label, dist.Normal(0,0.5))
    bA = norm('bA')
    bC = norm('bC')
    bI = norm('bI')
    bIA = norm('bIA')
    bIC = norm('bIC')
    bAge = norm('bAge')
    bEdu = norm('bEdu')
    edu_coefs = numpyro.sample('edu_coefs', dist.Dirichlet(jnp.repeat(2,7)))
    bMale = norm('bMale')
    # likelihood
    delta_j = jnp.pad(edu_coefs, (1, 0))
    delta_E = jnp.sum(jnp.where(jnp.arange(8) <= edu_cat[..., None], delta_j, 0), -1)
    BI = bI + bIA*action + bIC*contact
    phi = bA*action + BI*intention + bC*contact + bAge*age_std + bEdu*delta_E + bMale*male
    response_hat = numpyro.sample('response_hat', dist.OrderedLogistic(phi, cutpoints), obs=response-1)
kernel = NUTS(model)
mcmc = MCMC(kernel, num_warmup=500, num_samples=500, num_chains=4, chain_method='sequential')
dat = {k:v.to_numpy() for k,v in d[['action', 'contact', 'intention', 'age_std', 'edu_cat', 'male', 'response']].items()}
mcmc.run(rng_key, **dat)
mcmc.print_summary()<jupyter_output>sample: 100%|██████████| 1000/1000 [01:41<00:00,  9.86it/s, 63 steps of size 8.29e-02. acc. prob=0.96]
sample: 100%|██████████| 1000/1000 [01:11<00:00, 14.00it/s, 31 steps of size 8.61e-02. acc. prob=0.95]
sample: 100%|██████████| 1000/1000 [01:08<00:00, 14.54it/s, 31 steps of size 8.76e-02. acc. prob=0.95]
sample: 100%|██████████| 1000/1000 [01:18<00:00, 12.67it/s, 63 steps of size 8.40e-02. acc. prob=0.95]
{'fraction_non_alphanumeric': 0.13723061679464948, 'fraction_numerical': 0.04570225414912064, 'mean_word_length': 3.8237753882915175, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 15, 'https://': 7, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '533335', 'license': 'MIT', 'n_tokens_mistral': 3039, 'n_tokens_neox': 2681, 'n_words': 839, 'path': '_statrethink_numpyro_w07.ipynb', 'repo_name': 'robom/statrethink-course-numpyro-2019', 'size': '8074'}