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

# R packages 

```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(ggpubr)
library(cmdstanr)
library(bayesplot)
library(glmmTMB)
```

# Data 


```{r}
dat <- readxl::read_xlsx("../../data_confid/Behavior_Counts.xlsx") |> 
  mutate(calfID = as.numeric(as.factor(calfID)),
         timepoint_n = timepoint+2,
         timepoint_f  = as.numeric(as.factor(timepoint+2)))

# viz
dat |> 
  ggplot(aes(timepoint, behavior_count))+
  geom_point(aes(group = breed, 
                 fill = breed),
             shape =21, size= 2.5)+
  facet_wrap(~trt)+
  scico::scale_fill_scico_d()+
  theme_pubclean()
```

```{r}
m2 <- glmmTMB(behavior_count ~ timepoint * trt * breed +
                toep(timepoint + 0 | calfID),
              ziformula = ~1,
              family = poisson(link = "log"),
              REML  =T,
              data = dat |> mutate(timepoint_f=as.factor(timepoint)))
```



```{r}
model <- cmdstan_model("desexp.stan")

X <- model.matrix(behavior_count ~ 0 + timepoint * trt * breed, data = dat)

stan_data <- list(
  N = nrow(dat), 
  behavior_count = dat$behavior_count, 
  K_mu = ncol(X),
  X_mu = X, 
  J_calf = n_distinct(dat$calfID), 
  calfID = dat$calfID, 
  T_max = n_distinct(dat$timepoint_f), 
  timepoint_f = dat$timepoint_f)

fitted <- 
  model$sample(
    data = stan_data,
    seed = 42, 
    chains = 3, 
    parallel_chains = 3)

```






```{r warning=FALSE}
model2 <- cmdstan_model("desexp_v02.stan")

X <- model.matrix(behavior_count ~ 0 + factor(timepoint_f) * trt * breed, data = dat)

stan_data <- list(
  N = nrow(dat), 
  behavior_count = dat$behavior_count, 
  K_mu = ncol(X),
  X_mu = X, 
  J_calf = n_distinct(dat$calfID), 
  calfID = dat$calfID, 
  T_max = n_distinct(dat$timepoint_f), 
  timepoint_f = dat$timepoint_f, 
  times = unique(dat$timepoint_n))

fitted_v02 <- 
  model2$sample(
    data = stan_data,
    seed = 42, 
    chains = 3, 
    parallel_chains = 3, 
    iter_warmup = 1500, 
    iter_sampling = 2500)

```



```{r}
bayesplot::mcmc_trace(fitted$draws(c("alpha_zi", "beta_mu")))
bayesplot::mcmc_trace(fitted_v02$draws(c("phi", "beta_mu[1]", "beta_mu[10]")))
```

```{r}
hist(rnorm(1000, 0, sqrt(5)))
hist(rbeta(1000, 1, 1))
```

```{r}
colnames(X)
```

```{r}
bayesplot::mcmc_hist(fitted_v02$draws(c("phi", "beta_mu")))
```

```{r}
hist(rnorm(1000, 0, sqrt(.5)))

bayesplot::mcmc_hist(fitted_v02$draws(c("rho_toep")))
```











