NASA_chi_square_test.Rmd

# 航天事业的预算有党派门户之见 {#NASA}


```{r}
library(tidyverse)
library(tidybayes)
library(bayesplot)
library(rstan)
library(loo)
library(broom.mixed)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
theme_set(bayesplot::theme_default())
```


## 导入数据
```{r}
gss <- read_rds("./rawdata/gss.rds")
```


```{r}
gss %>%
  select(party, NASA) %>%
  ggplot(aes(x = party, y = NASA, color = NASA)) +
  geom_jitter()
```



```{r}
gss %>%
  select(party, NASA) %>%
  ggplot(aes(x = party, fill = NASA)) +
  geom_bar()
```



```{r}
chisq.test(gss$party, gss$NASA)
```


## using the categorical distribution



```{r, warning=FALSE, message=FALSE}
stan_program <- "
data {
  int<lower=1> N;        
  int<lower=2> C;            
  int<lower=1, upper=C> y[N]; 
}
parameters {
  simplex[C] theta;
}
model {
  target += dirichlet_lpdf(theta | rep_vector(2, C));
  for(n in 1:N) {
    target += categorical_lpmf(y[n] | theta);
  }
}

"


stan_data <- gss %>% 
  tidybayes::compose_data(
    N     = nrow(gss), 
    C     = n_distinct(NASA),
    y     = gss$NASA
)

stan_fit0 <- stan(model_code = stan_program, data = stan_data,
                  iter = 1000, chains = 4)
```

```{r}
stan_fit0
```

## using the categorical distribution and grouping

```{r, warning=FALSE, message=FALSE}
stan_program <- "
data {
  int<lower=1> N;        
  int<lower=2> C;            
  int<lower=2> J;            
  int<lower=1, upper=C> y[N]; 
  int<lower=1, upper=J> g[N]; 
}
parameters {
  simplex[C] theta[J];
}
model {
  for(j in 1:J) {
    target += dirichlet_lpdf(theta[j] | rep_vector(2, C));
  }
  
  for(n in 1:N) {
    target += categorical_lpmf(y[g[n]] | theta[g[n]]);
  }
}
generated quantities {
  matrix[C, J-1] diff;
  for(i in 1:C) {
    for(j in 1:J-1) {
      diff[i, j] = theta[i, j] - theta[i, j+1];
    }
  }
}
"


stan_data <- gss %>% 
  tidybayes::compose_data(
    N     = nrow(gss), 
    C     = n_distinct(NASA),
    J     = n_distinct(party),
    g     = gss$party, 
    y     = gss$NASA
)

stan_fit <- stan(model_code = stan_program, data = stan_data,iter = 1000, chains = 4)
```


```{r}
stan_fit
```


```{r}
stan_fit %>% 
  tidybayes::gather_samples(diff[i,j]) %>% # catogroy  - group1- group2
  group_by(i, j) %>%
  mutate(jj = if_else(j == 1, "group1_minus_group2", "group2_minus_group3")) %>%
  mutate(diff = str_c("cate", i, jj, sep = "_")) %>%
  ggplot(aes(x = estimate, y = diff)) +
  tidybayes::geom_halfeyeh() +
  labs(y = NULL)
```


## using multimonial

```{r}
d <- gss %>%
  select(NASA, party) %>%
  count(party, NASA) 
d
```





```{r}
d <- gss %>%
  select(NASA, party) %>%
  count(NASA) %>%
  mutate(sum = sum(n))
d
```


- 与categorical不同的是，multinomial 需要弄成表格的形式


https://vasishth.github.io/bayescogsci/book/modeling-multiple-categorical-responses.html

```{r, warning=FALSE, message=FALSE}
stan_program <- "
data {
  int<lower=2> C;  
  int<lower=0> y[C];  
}
parameters {
  simplex[C] theta;
}
model {
  target += dirichlet_lpdf(theta | rep_vector(2, C));
  y ~ multinomial(theta); //int[ ] ~ multinomial(vector) 

}

"


stan_data <- list(
    #N_trials = unique(d$sum),
    C        = length(unique(d$NASA)),
    y        = d$n
)

stan_eq <- stan(model_code = stan_program, data = stan_data,iter = 1000, chains = 4)
```




## revision edition


```{r}
d2 <- gss %>%
  select(NASA, party) %>%
  count(party,NASA) 

d2
```

```{r, warning=FALSE, message=FALSE}
stan_program <- "
data {
  int<lower=2> C;  
  int<lower=2> J;  
  vector[C] y[J];  
}
parameters {
  simplex[C] theta[J];
}
model {
  target += dirichlet_lpdf(theta | rep_vector(2, C));
  for(j in 1:J) {
    target += multinomial_lpmf(y[j] | theta[j]);
  } 

}

"


stan_data <- list(
    #N_trials = unique(d$sum),
    C        = length(unique(d$NASA)),
    y        = d$n
)

stan_eq <- stan(model_code = stan_program, data = stan_data,iter = 1000, chains = 4)
```