---
title: "Analysis of the Greek data"
author: "Michel RAYMOND"
output:
  html_document:
    toc: yes
    toc_float: yes
    number_section: no
    df_print: paged
    dev: png
  pdf_document:
    toc: yes
    number_section: yes
    keep_tex: yes
    highlight: "tango"
    df_print: kable
    dev: png
---
<style type="text/css">

body{ /* Normal  */
      font-size: 12px;
  }
td {  /* Table  */
  font-size: 14px;
}
h1.title {
  font-size: 38px;
  color: DarkRed;
}
h1 { /* Header 1 */
  font-size: 28px;
  color: DarkBlue;
}
h2 { /* Header 2 */
    font-size: 22px;
  color: DarkBlue;
}
h3 { /* Header 3 */
  font-size: 18px;
  font-family: "Times New Roman", Times, serif;
  color: DarkBlue;}
h4 { /* Header 4 */
  font-size: 16px;
  font-family: "Times New Roman", Times, serif;
  color: Blue;
}
code.r{ /* Code block */
    font-size: 12px;
}
pre { /* Code block - determines code spacing between lines */
    font-size: 14px;
}
</style>
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```

```{r library,          include = FALSE,message = FALSE}
library(car)
library(nimble)
library(basicMCMCplots)
library(MCMCglmm)
library(kableExtra)
library(knitr) 
```
```{r nimble_functions, include = FALSE,message = FALSE}            
        sampler_RJ_indicator2 <- nimbleFunction(
          name = 'sampler_RJ_indicator2',
          contains = sampler_BASE,
          setup = function(model, mvSaved, target, control) {
            coefNode1 <- control$coefNode1
            coefNode2 <- control$coefNode2
            propScale1 <- control$propScale1
            propScale2 <- control$propScale2
            calcNodes <- model$getDependencies(c(coefNode1, coefNode2, target))
            calcNodesReduced <- model$getDependencies(target)
          },
          run = function() {
            currentIndicator <- model[[target]]
            if(currentIndicator == 0) {
              currentLP <- model$getLogProb(calcNodesReduced)
              propCoef1 <- rnorm(1, 0, propScale1)
              propCoef2 <- rnorm(1, 0, propScale2)
              forwardPropLP <- dnorm(propCoef1, 0, propScale1, log = TRUE) + dnorm(propCoef2, 0, propScale2, log = TRUE)
              model[[target]] <<- 1
              model[[coefNode1]] <<- propCoef1
              model[[coefNode2]] <<- propCoef2
              proposalLP <- model$calculate(calcNodes)
              logMHR <- proposalLP - currentLP - forwardPropLP
            } else {
              currentLP <- model$getLogProb(calcNodes)
              currentCoef1 <- model[[coefNode1]]
              currentCoef2 <- model[[coefNode2]]
              reversePropLP <- dnorm(currentCoef1, 0, propScale1, log = TRUE) + dnorm(currentCoef2, 0, propScale2, log = TRUE)
              model[[target]] <<- 0
              model[[coefNode1]] <<- 0
              model[[coefNode2]] <<- 0
              model$calculate(calcNodes)
              proposalLP <- model$getLogProb(calcNodesReduced)
              logMHR <- proposalLP - currentLP + reversePropLP
            }
            accept <- decide(logMHR)
            if(accept) { copy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
            } else     { copy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE) }
          },
          methods = list(reset = function() { })
        )
        
        sampler_RJ_toggle2 <- nimbleFunction(
          name = 'sampler_RJ_toggle2',
          contains = sampler_BASE,
          setup = function(model, mvSaved, target, control) {
            indicatorNode <- control$indicatorNode
            samplerConf <- control$samplerConf
            samplerFL <- nimbleFunctionList(sampler_BASE)
            samplerFL[[1]] <- samplerConf$buildSampler(model = model, mvSaved = mvSaved)
          },
          run = function() { if(model[[indicatorNode]] != 0)   samplerFL[[1]]$run() },
          methods = list(
            reset = function() { samplerFL[[1]]$reset() }
          )
        )
        sampler_RJ_indicator1 <- nimbleFunction(   
          name = 'sampler_RJ_indicator1',
          contains = sampler_BASE,
          setup = function(model, mvSaved, target, control) {
            coefNode1 <- control$coefNode1
            propScale1 <- control$propScale1
            calcNodes <- model$getDependencies(c(coefNode1, target))  
            calcNodesReduced <- model$getDependencies(target)
          },
          run = function() {
            currentIndicator <- model[[target]]
            if(currentIndicator == 0) {
              currentLP <- model$getLogProb(calcNodesReduced)
              propCoef1 <- rnorm(1, 0, propScale1)
              forwardPropLP <- dnorm(propCoef1, 0, propScale1, log = TRUE) 
              model[[target]] <<- 1
              model[[coefNode1]] <<- propCoef1
              proposalLP <- model$calculate(calcNodes)
              logMHR <- proposalLP - currentLP - forwardPropLP
            } else {
              currentLP <- model$getLogProb(calcNodes)
              currentCoef1 <- model[[coefNode1]]
              reversePropLP <- dnorm(currentCoef1, 0, propScale1, log = TRUE) 
              model[[target]] <<- 0
              model[[coefNode1]] <<- 0
              model$calculate(calcNodes)
              proposalLP <- model$getLogProb(calcNodesReduced)
              logMHR <- proposalLP - currentLP + reversePropLP
            }
            accept <- decide(logMHR)
            if(accept) { copy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
            } else     { copy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE) }
          },
          methods = list(reset = function() { })
        )
```

```{r read_file}
            file_name<-"Greek_data.csv"
            df<-read.csv(file_name)
            df$yob_s<-as.numeric(scale(df$yob, center = TRUE, scale=T))
```
File name:  `r file_name` 

## Descriptive statistics

```{r , echo=FALSE}
            N       <- nrow(df)
            Nsib    <- df$Nsib
            sexPref <- df$sexPref
            os      <- df$os
            ob      <- df$ob
            yob     <- df$yob
            age     <- df$age
ds<-data.frame(item=NA,hetero=NA,homo=NA,All=NA)
# N
  N_hetero<-length(which(df$sexPref==0))
  N_homo  <-length(which(df$sexPref==1))
  l1<-c('N',N_hetero,N_homo,N)
  l2<-c(' %',round((N_hetero/N)*100,1),round((N_homo/N)*100,1),100)
  ds<-rbind(ds,l1,l2)
#Age
  age_hetero<-round(as.numeric(tapply(age, sexPref, mean)[1]),1)
  age_homo  <-round(as.numeric(tapply(age, sexPref, mean)[2]),1)
  l3<-c('Age (mean)',age_hetero,age_homo, round(mean(age),1))
  ds<-rbind(ds,l3)
  age_h_sd<-round(as.numeric(tapply(age, sexPref, sd)[1]),1)
  age_H_sd<-round(as.numeric(tapply(age, sexPref, sd)[2]),1)
  l4<-c('Age   (sd)',age_h_sd,age_H_sd,round(sd(age),1))
   ds<-rbind(ds,l4)
#Yob
  yob_hetero<-round(as.numeric(tapply(yob, sexPref, mean)[1]),1)
  yob_homo  <-round(as.numeric(tapply(yob, sexPref, mean)[2]),1)
  l31<-c('Yob (mean)',yob_hetero,yob_homo, round(mean(yob),1))
  ds<-rbind(ds,l31)
# Nsib (sum)
   f_hetero<-round(sum(Nsib[which(sexPref==0)]),2)
   f_homo  <-round(sum(Nsib[which(sexPref==1)]),2)
   f_all   <-round(sum(Nsib),2)
   l5<-c('Nsib (sum)',f_hetero,f_homo ,f_all)  
    ds<-rbind(ds,l5)
# Nsib (mean)
   f_hetero<-round(mean(Nsib[which(sexPref==0)]),2)
   f_homo  <-round(mean(Nsib[which(sexPref==1)]),2)
   f_all   <-round(mean(Nsib),2)
   l5<-c('Nsib (mean)',f_hetero,f_homo ,f_all)  
    ds<-rbind(ds,l5)
# Older brothers (mean)   
   ob_hetero<-round(mean(ob[which(sexPref==0)]),2)
   ob_homo  <-round(mean(ob[which(sexPref==1)]),2)
   ob_all   <-round(mean(ob),2)
   l6<-c('OB (mean)',ob_hetero,ob_homo ,ob_all) ; ds<-rbind(ds,l6) 

# Older sisters
   os_hetero<-round(mean(os[which(sexPref==0)]),2)
   os_homo  <-round(mean(os[which(sexPref==1)]),2)
   os_all   <-round(mean(os),2)
   l7<-c('OS (mean)',os_hetero,os_homo ,os_all) ; ds<-rbind(ds,l7) 

  ds<-ds[-1,]
  names(ds)[1]<-" "
  kbl(ds) %>%
  kable_classic(full_width = T, html_font = "Cambria", font_size = 20) 
```

## Mothers: Sexual orientation and fertility (controlled for birth rank)

### Birth rank as a qualitative variable

```{r test_mother_quali}
              sexPref <- df$sexPref
              ob      <- df$ob
              os      <- df$os
              Nsib    <- df$Nsib
              yob_s   <- df$yob_s
              sexPref<-as.factor(sexPref)
              sexPref<-relevel(sexPref, ref='0')
              rang<-ob+1
              rang<-as.factor(rang)
              rang<-relevel(rang, ref=1)
              sexPref<-as.factor(sexPref)
              summary(m1<-glm(Nsib~sexPref+rang+yob_s, family="poisson")        )
              table<-Anova(m1)
              table %>%
              kbl(caption = "Anova type II") %>%
              kable_classic(full_width = F, html_font = "Cambria")   
```

### Birth rank as a quantitative variable

```{r test_mother_quanti}
              rang<-as.character(rang)
              rang<-as.numeric(rang)
              summary(m1_<-glm(Nsib~sexPref+rang+yob_s, family="poisson")        )
              table<-Anova(m1_)
              table %>%
              kbl(caption = "Anova type II") %>%
              kable_classic(full_width = F, html_font = "Cambria")  
```


## Analysis
```{r presentation_models}

    model<-paste('model',c(2,4,5,6),sep='')
     OBE<-c('  x','  x', '  x' ,'{x}' )
      AE<-c(''  ,'  x', '{x}', '  x' )
Comments<-c('OBE estimates','OBE & AE estimates','AE support (in presence of OBE)' ,'OBE support (in presence of AE)')
fo<-data.frame(Name=model,OBE=OBE,AE=AE,Purpose=Comments)
                fo %>%
                kbl(caption = "The various models considered") %>%
                kable_classic(full_width = F, html_font = "Cambria")   #library(kableExtra)
```

## Older sibs effect alone

### WAIC comparison 

#### Data restriction 

To avoid the effects of small sample sizes for the total number of older brothers (OB) or sisters (OS) in categories poorly sampled, data are restricted to categories (total number of OB or OS) including at least 10 individuals.
 
```{r data_restriction}
print('Number of sampled individuals with 0, 1, 2, 3, 4, or OB');t1<-table(df$sexPref,df$ob)
rownames(t1) = c("hetero", "homo") 
t1
print(paste('For OB: ',length(which(df$ob>2)),' individuals suppressed with ob>2'))

print('Number of sampled individuals with 0, 1, 2, 3, or 9 OS');t2<-table(df$sexPref,df$os)
rownames(t2) = c("hetero", "homo") 
t2
#suppression des individus avec ob>2
print(paste('For OS: ',length(which(df$os>2)),' individuals suppressed with os>2'))

df_ob<-df[-which(df$ob>2),]
df_os<-df[-which(df$os>2),]
```

#### Chain parameters

```{r WAIC_init, echo=TRUE, eval=TRUE}
    nrep <- 10
    nchains <- 5
    niter <-   50000
    nburnin <- 20000
```
```{r VariousOBEfunctionsOBE,  include = FALSE,message = TRUE}
            N       <- nrow(df_ob)
            sexPref <- df_ob$sexPref
            ob      <- df_ob$ob
modelInfoList <- list(
    ##  ilogit(a*B + b)*c+p0  ####
    `f1` = list(
         code = nimbleCode({
                  a ~ dunif(-100, 100)
                  b ~ dunif(-50, 50)
                  c ~ dunif(0,1)
                  p0 ~ dunif(0, 1)
                  for(i in 1:N) {
                    p[i] <- min(ilogit(a*ob[i] + b)*c+p0 , 0.999999)
                    sexPref[i] ~ dbern(p[i])
                  }
                }),
        constants = list(N=N, ob=ob),
        data = list(sexPref=sexPref),
        inits = list(a=0.1,b=-1,p0=0.5,c=0.5)
    ),
    ##  ilogit(a*B/sib + b)*(1-p0)+p0  ####
    `f2` = list(
         code = nimbleCode({
                  a ~ dunif(-100, 100)
                  b ~ dunif(-50, 50)
                  p0 ~ dunif(0, 1)
                  for(i in 1:N) {
                    p[i] <- ilogit(a*ob[i] + b)*(1-p0)+p0
                    sexPref[i] ~ dbern(p[i])
                  }
                }),
        constants = list(N=N, ob=ob),
        data = list(sexPref=sexPref),
        inits = list(a=0.1,b=-1,p0=0.5)
    ),
    ##  ilogit(a*B/sib + logit(p0))  ####
    `f3` = list(
         code = nimbleCode({
                  a ~ dunif(-100, 100)
                  p0 ~ dunif(0, 1)
                  b<-logit(p0)
                  for(i in 1:N) {
                    p[i] <- ilogit(a*ob[i] + b)
                    sexPref[i] ~ dbern(p[i])
                  }
                }),
        constants = list(N=N, ob=ob),
        data = list(sexPref=sexPref),
        inits = list(a=0.1,p0=0.5)
    ),   
    ##  a*B^g   ####
    `f4` = list(
         code = nimbleCode({
          p0 ~ dunif(0, 1)
          a ~ dunif(-5, 5)
          gamma~ dunif(0, 20) 
          for(i in 1:N) {
             p[i] <- min(p0 +  a*(ob[i])^gamma, 0.999999)
             sexPref[i] ~ dbern(p[i])
            }
        }),
        constants = list(N=N,ob=ob),
        data = list(sexPref=sexPref),
        inits = list(p0=0.5, a=0.2,gamma=0.1)
    ),
    ##  line(B)  ####
    `f5` = list(
         code = nimbleCode({
                  a ~ dunif(-100, 100)
                  p0 ~ dunif(0,1)
                  for(i in 1:N) {
                    p[i] <- min(a*ob[i]+p0 , 0.999999)
                    sexPref[i] ~ dbern(p[i])
                  }
                }),
        constants = list(N=N, ob=ob),
        data = list(sexPref=sexPref),
        inits = list(a=0.1,p0=0.5)
    ),   
    ##  p0, p0+a1 + line(B)  ####
    `f6` = list(
         code = nimbleCode({
            a ~ dunif(-100, 100)
            a1 ~ dunif(0,1)
            p0 ~ dunif(0,1)
            for(i in 1:N) {
              p[i]       <- min(p0 + step(ob[i]-1)*(a1-a+a*ob[i]),0.99999) 
              sexPref[i] ~ dbern(p[i])
            }
        }),
        constants = list(N=N, ob=ob),
        data = list(sexPref=sexPref),
        inits = list(a=0.1,p0=0.5,a1=0.1)
    ),   
    ##  p0, p0+a1   ####
    `f7` = list(
         code = nimbleCode({
            a1 ~ dunif(0,1)
            p0 ~ dunif(0,1)
            for(i in 1:N) {
              p[i]       <- min(p0 + step(ob[i]-1)*a1,0.99999) 
              sexPref[i] ~ dbern(p[i])
            }
        }),
        constants = list(N=N, ob=ob),
        data = list(sexPref=sexPref),
        inits = list(p0=0.5,a1=0.1)
      ),
    ##  polynomiale ####
    `f8` = list(
         code = nimbleCode({
            a2 ~ dunif(-100, 100)
            a1 ~ dunif(-100, 100)
            p0 ~ dunif(0,1)
            for(i in 1:N) {
              p[i]       <- min(p0 + a1*ob[i]+a2*ob[i]^2,0.99999) 
              sexPref[i] ~ dbern(p[i])
            }
        }),
        constants = list(N=N, ob=ob),
        data = list(sexPref=sexPref),
        inits = list(a1=0.1,a2=0.1,p0=0.5)
    )
)
    ### run ####

runModel2 <- function(mi, niter, nburnin, nchains, nrep) {
    Rmodel <- nimbleModel(mi$code, mi$constants, mi$data, mi$inits)
    lp <- Rmodel$calculate()
    if(is.na(lp) || is.nan(lp)) stop('model not initialized')
    conf <- configureMCMC(Rmodel, enableWAIC = TRUE)
    Rmcmc <- buildMCMC(conf)
    compiledList <- compileNimble(list(model=Rmodel, mcmc=Rmcmc))
    Cmodel <- compiledList$model
    Cmcmc <- compiledList$mcmc
    waics <- numeric(nrep)
    for(i in 1:nrep) {
        Cmodel$setInits(mi$inits)
        waics[i] <- runMCMC(Cmcmc,
                            niter=niter, nburnin=nburnin, nchains=nchains,
                            progressBar = FALSE, samples = FALSE,
                            summary = FALSE, WAIC = TRUE)
    }
    return(waics)
}
    set.seed(round(runif(1, min=0, max=1000), digits=0))
    waicList2 <- lapply(modelInfoList, runModel2, niter, nburnin, nchains, nrep)
    dfobe <- do.call('rbind',
                   lapply(1:length(waicList2),
                          function(i) data.frame(Model = names(waicList2)[i],
                                                 WAIC = waicList2[[i]])
                          )
                   )
```
```{r VariousOBEfunctionsOSE,  include = FALSE,message = TRUE}
            N       <- nrow(df_os)
            sexPref <- df_os$sexPref
            os      <- df_os$os
modelInfoList <- list(
    ##  ilogit(a*S+ b)*c+p0   ####
    `f2` = list(
         code = nimbleCode({
                  a ~ dunif(-100, 100)
                  b ~ dunif(-50, 50)
                  c ~ dunif(0,1)
                  p0 ~ dunif(0, 1)
                  for(i in 1:N) {
                    p[i] <- min(ilogit(a*os[i] + b)*c+p0 , 0.999999)
                    sexPref[i] ~ dbern(p[i])
                  }
                }),
        constants = list(N=N,os=os),
        data = list(sexPref=sexPref),
        inits = list(a=0.1,b=-1,p0=0.5,c=0.5)
    ),
    ##  ilogit(a*S/sib + b)*(1-p0)+p0  ####
    `f5` = list(
         code = nimbleCode({
                  a ~ dunif(-100, 100)
                  b ~ dunif(-50, 50)
                  p0 ~ dunif(0, 1)
                  for(i in 1:N) {
                    p[i] <- ilogit(a*os[i] + b)*(1-p0)+p0
                    sexPref[i] ~ dbern(p[i])
                  }
                }),
        constants = list(N=N, os=os),
        data = list(sexPref=sexPref),
        inits = list(a=0.1,b=-1,p0=0.5)
    ),
    `f8` = list(
         code = nimbleCode({
                  a ~ dunif(-100, 100)
                  p0 ~ dunif(0, 1)
                  b<-logit(p0)
                  for(i in 1:N) {
                    p[i] <- ilogit(a*os[i] + b)
                    sexPref[i] ~ dbern(p[i])
                  }
                }),
        constants = list(N=N, os=os),
        data = list(sexPref=sexPref),
        inits = list(a=0.1,p0=0.5)
    ),    
    ##  a*S^g  ####
    `f14` = list(
         code = nimbleCode({
          p0 ~ dunif(0, 1)
          a ~ dunif(-5, 5)
          gamma~ dunif(0, 20) 
          for(i in 1:N) {
             p[i] <- min(p0 +  a*(os[i])^gamma, 0.999999)
             sexPref[i] ~ dbern(p[i])
            }
        }),
        constants = list(N=N,os=os),
        data = list(sexPref=sexPref),
        inits = list(p0=0.5, a=0.2,gamma=0.1)
    ),
    ##  line(S)  ####
    `f20` = list(
         code = nimbleCode({
                  a ~ dunif(-100, 100)
                  p0 ~ dunif(0,1)
                   for(i in 1:N) {
                    p[i] <- min(a*os[i]+p0 , 0.999999)
                    sexPref[i] ~ dbern(p[i])
                  }
                }),
        constants = list(N=N, os=os),
        data = list(sexPref=sexPref),
        inits = list(a=0.1,p0=0.5)
    ),
    ## p0, p0+a1 + line(S)  ####
    `f23` = list(
         code = nimbleCode({
             a ~ dunif(-100, 100)
            a1 ~ dunif(0,1)
            p0 ~ dunif(0,1)
            for(i in 1:N) {
              p[i]       <- min(p0 + step(os[i]-1)*(a1-a+a*os[i]),0.99999) 
              sexPref[i] ~ dbern(p[i])
                  }
                }),
        constants = list(N=N, os=os),
        data = list(sexPref=sexPref),
        inits = list(a=0.1,p0=0.5,a1=0.1)
    ),
    ## p0, line(S)  ####
    `f26` = list(
         code = nimbleCode({
            a1 ~ dunif(0,1)
            p0 ~ dunif(0,1)
            for(i in 1:N) {
              p[i]       <- min(p0 + step(os[i]-1)*a1,0.99999) 
              sexPref[i] ~ dbern(p[i])
                  }
                }),
        constants = list(N=N, os=os),
        data = list(sexPref=sexPref),
        inits = list(p0=0.5,a1=0.1)
    ),
    ## polynomiale(S)  ####
    `f29` = list(
         code = nimbleCode({
            a2 ~ dunif(-100, 100)
            a1 ~ dunif(-100, 100)
            p0 ~ dunif(0,1)
            for(i in 1:N) {
              p[i]       <- min(p0 + a1*os[i]+a2*os[i]^2,0.99999) 
              sexPref[i] ~ dbern(p[i])
                  }
                }),
        constants = list(N=N, os=os),
        data = list(sexPref=sexPref),
        inits = list(a1=0.1,a2=0.1,p0=0.5)
    )
)
    ### run ####
runModel2 <- function(mi, niter, nburnin, nchains, nrep) {
    Rmodel <- nimbleModel(mi$code, mi$constants, mi$data, mi$inits)
    lp <- Rmodel$calculate()
    if(is.na(lp) || is.nan(lp)) stop('model not initialized')
    conf <- configureMCMC(Rmodel, enableWAIC = TRUE)
    Rmcmc <- buildMCMC(conf)
    compiledList <- compileNimble(list(model=Rmodel, mcmc=Rmcmc))
    Cmodel <- compiledList$model
    Cmcmc <- compiledList$mcmc
    waics <- numeric(nrep)
    for(i in 1:nrep) {
        Cmodel$setInits(mi$inits)
        waics[i] <- runMCMC(Cmcmc,
                            niter=niter, nburnin=nburnin, nchains=nchains,
                            progressBar = FALSE, samples = FALSE,
                            summary = FALSE, WAIC = TRUE)
    }
    return(waics)
}
    set.seed(round(runif(1, min=0, max=1000), digits=0))
    waicList2 <- lapply(modelInfoList, runModel2, niter, nburnin, nchains, nrep)
    dfose <- do.call('rbind',
                   lapply(1:length(waicList2),
                          function(i) data.frame(Model = names(waicList2)[i],
                                                 WAIC = waicList2[[i]])
                          )
                   )
    
```
```{r tableWAIC,             include=TRUE, eval=TRUE}
model<-c('f1','f2',"f3","f4","f5","f6","f7","f8")
#OB
  dfobe$Model<-factor(dfobe$Model)
  ob_waic<-tapply(dfobe$WAIC,dfobe$Model,'mean')
  
fo<-data.frame(Function=model,WAIC=ob_waic)
   row.names(fo)<-NULL
   fo %>%
   kbl(caption = "WAIC for OB") %>%
   kable_classic(full_width = F, html_font = "Cambria") 
#OS
  dfose$Model<-factor(dfose$Model)
  os_waic<-tapply(dfose$WAIC,dfose$Model,'mean')

  fo<-data.frame(Function=model,WAIC=os_waic)
  row.names(fo)<-NULL
   fo %>%
   kbl(caption = "WAIC for OS") %>%
   kable_classic(full_width = F, html_font = "Cambria") 
```

## Estimates: AE & OBE (model 4)

```{r}
  N       <- nrow(df_ob)
  sexPref <- df_ob$sexPref
  ob      <- df_ob$ob
  yob     <- df_ob$yob_s 
```
```{r model4_f25 ,echo=TRUE       }
 niter <-  100000
 nburnin <- 10000
 code <- nimbleCode({
           a1 ~ dunif(0,1)
      phi_ea ~ dunif(0, 1)
         phi ~ dunif(0, 1)
          p0 ~ dunif(0, 1)
       alpha ~ dunif(0 , 2)
        beta ~ dunif(0, 10)
        c[1] ~ dunif(-2, 10)
        c[2] ~ dunif(-10, 10)
  beta_hyper ~ dunif(0, 10)
  for(i in 1:N) {
       lambda[i] <- exp(c[1] + c[2]*yob[i])*(1+beta*x[i])*(1+beta_hyper*x_hy[i])
            p1[i]<-min((1+alpha*x[i])*p0,0.999999)
             p[i]<- min(p1[i] + step(ob[i]-1)*a1,0.99999) 
         Nsib[i] ~ dpois(lambda[i])       
      sexPref[i] ~ dbern(p[i])
            x[i] ~ dbern(phi_ea)
         x_hy[i] ~ dbern(phi)
  }
})
constants = list(N=N,yob=yob,ob=ob)
data = list(Nsib=Nsib,sexPref=sexPref)
inits = list(a1=0.2, c=c(1,0.1), x = rbinom(N, size=1, prob=0.5),
              x_hy = rbinom(N, size=1, prob=0.5), beta_hyper = 1,
              phi_ea = .5, phi=0.5,alpha = .5, beta = .5,p0=0.3)
```
```{r model4_f25_suite, eval=TRUE}
                      Rmodel <- nimbleModel(code, constants, data, inits)
                      nimbleOptions(showCompilerOutput = TRUE) 
                      lp <- Rmodel$calculate()
                      if(is.na(lp) || is.nan(lp)) stop('model not initialized')
                      Rmodel$initializeInfo()
                      conf <- configureMCMC(Rmodel)
                      Rmcmc <- buildMCMC(conf)
                      Cmodel <- compileNimble(Rmodel)
                      Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
                      set.seed(round(runif(1, min=0, max=1000), digits=0))
                      samples <- runMCMC(Cmcmc, niter, nburnin)
```
```{r model4_f25_figures ,eval=TRUE}
               table<-samplesSummary(samples)
                mode<-NULL
                for (j in 1:dim(table)[1]) {
                  mode[j]  <-mean(density(samples[,j])$x[which(density(samples[,j])$y==max(density(samples[,j])$y))])# mode of the posterior distribution
                }
                table_out<-round(cbind(table,mode),3)
                table_out %>%
                kbl(caption = "Estimates model 4 with function 7 (OB)") %>%
                kable_classic(full_width = F, html_font = "Cambria")   
```

## Support:  AE in presence of an OBE (model 5)
 

```{r model5_6_OBE_init, echo=FALSE}  
    N       <- nrow(df_ob)
    sexPref <- df_ob$sexPref
    Nsib    <- df_ob$Nsib
    ob      <- df_ob$ob
    yob     <- df_ob$yob_s 
```

### OBE function 7

```{r model5_obe_f7, echo = TRUE}
niter <- 500000
nburnin <- 10
code <- nimbleCode({
           a1 ~ dunif(0,1)
      phi_ea ~ dunif(0, 1)
         phi ~ dunif(0, 1)
          p0 ~ dunif(0, 1)
       alpha ~ dunif(0 , 2)
        beta ~ dunif(0, 10)
        c[1] ~ dunif(-2, 10)
        c[2] ~ dunif(-10, 10)
        ant ~ dbern(0.5)
  beta_hyper ~ dunif(0, 10)
  for(i in 1:N) {
       lambda[i] <- exp(c[1] + c[2]*yob[i])*(1+beta*x[i]*ant)*(1+beta_hyper*x_hy[i])
       p1[i]<-min((1+alpha*x[i]*ant)*p0,0.999999)
       p[i]<- min(p1[i] + step(ob[i]-1)*a1,0.99999) 
         Nsib[i] ~ dpois(lambda[i])       
      sexPref[i] ~ dbern(p[i])
            x[i] ~ dbern(phi_ea)
         x_hy[i] ~ dbern(phi)
  }
})
constants = list(N=N,yob=yob,ob=ob)
data = list(Nsib=Nsib,sexPref=sexPref)
inits = list(a1=0.1, c=c(1,0.1), x = rbinom(N, size=1, prob=0.5),
              x_hy = rbinom(N, size=1, prob=0.5), beta_hyper = 1,ant=0,
              phi_ea = .5, phi=0.5,alpha = .5, beta = .5,p0=0.3)
```
```{r model5_obe_f7_suite, eval=TRUE}
        
        Rmodel <- nimbleModel(code, constants, data, inits)
        Rmodel$calculate()
        
        conf <- configureMCMC(Rmodel, enableWAIC = TRUE)
        conf$printSamplers(byType = TRUE)
        conf$removeSampler('ant')
        conf$printSamplers(byType = TRUE)
        conf$addSampler(target = 'ant', type = 'RJ_indicator2',
                        control =
                          list(coefNode1 = 'alpha', propScale1 = 0.1,    
                               coefNode2 = 'beta',  propScale2 = 0.5)) 
        conf$printSamplers(byType = TRUE)
        nodes <- c('alpha', 'beta')  
        print(nodes)
        for(node in nodes) {
          samplerConf <- conf$getSamplers(node)[[1]]
          conf$removeSamplers(node)
          conf$addSampler(target = node, type = 'RJ_toggle2',
                          control =
                            list(indicatorNode = 'ant',
                                 samplerConf = samplerConf))
        }
        conf$printSamplers(byType = TRUE)
        conf$printMonitors()
        Rmcmc <- buildMCMC(conf)
        Cmodel <- compileNimble(Rmodel)
        Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
        set.seed(round(runif(1, min=0, max=1000), digits=0))
        system.time(samples <- runMCMC(Cmcmc, niter, nburnin))
```
```{r model5_obe_f7_densityplot, fig.height = 2, fig.align = "center",eval=TRUE}
                samplesPlot(samples, 'ant', densityplot = FALSE)
                samples_temp <- samples[1000:(niter-nburnin),]
                table<-samplesSummary(samples_temp)
                round(table,3) %>%
                kbl(caption = "Estimates") %>%
                kable_classic(full_width = F, html_font = "Cambria")
```

## Support: OBE in presence of AE (model 6)

### OBE function 7

```{r model6_obe_f7,echo=TRUE}
 niter <- 500000
 nburnin <- 10
 code <- nimbleCode({
           a1 ~ dunif(0,1)
      phi_ea ~ dunif(0, 1)
         phi ~ dunif(0, 1)
          p0 ~ dunif(0, 1)
       alpha ~ dunif(0 , 2)
        beta ~ dunif(0, 10)
        c[1] ~ dunif(-2, 10)
        c[2] ~ dunif(-10, 10)
        obe ~ dbern(0.5)
  beta_hyper ~ dunif(0, 10)
  for(i in 1:N) {
       lambda[i] <- exp(c[1] + c[2]*yob[i])*(1+beta*x[i])*(1+beta_hyper*x_hy[i])
       p1[i]<-min((1+alpha*x[i])*p0,0.999999)
         p[i]<- min(p1[i] + step(ob[i]-1)*a1*obe,0.99999) 
         Nsib[i] ~ dpois(lambda[i])       
      sexPref[i] ~ dbern(p[i])
            x[i] ~ dbern(phi_ea)
         x_hy[i] ~ dbern(phi)
  }
})
constants = list(N=N,yob=yob,ob=ob)
data = list(Nsib=Nsib,sexPref=sexPref)
inits = list(a1=0.1, c=c(1,0.1), x = rbinom(N, size=1, prob=0.5),
              x_hy = rbinom(N, size=1, prob=0.5), beta_hyper = 1,obe=0,
              phi_ea = .5, phi=0.5,alpha = .5, beta = .5,p0=0.3)
```
```{r model6_obe_f7_suite, eval=TRUE}
        
        Rmodel <- nimbleModel(code, constants, data, inits)
        Rmodel$calculate()
        conf <- configureMCMC(Rmodel, enableWAIC = TRUE)
        conf$printSamplers(byType = TRUE)
        conf$removeSampler('obe')
        conf$printSamplers(byType = TRUE)
        conf$addSampler(target = 'obe', type = 'RJ_indicator1',   
                                control =
                                  list(coefNode1 = 'a1', propScale1 = 0.3))  
        conf$printSamplers(byType = TRUE)
        nodes <- c('a1')  
        print(nodes)
        
        for(node in nodes) {
          samplerConf <- conf$getSamplers(node)[[1]]
          conf$removeSamplers(node)
          conf$addSampler(target = node, type = 'RJ_toggle2',
                          control =
                            list(indicatorNode = 'obe',
                                 samplerConf = samplerConf))
        }
        conf$printSamplers(byType = TRUE)
        conf$printMonitors()
        Rmcmc <- buildMCMC(conf)
        Cmodel <- compileNimble(Rmodel)
        Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
        set.seed(round(runif(1, min=0, max=1000), digits=0))
        system.time(samples <- runMCMC(Cmcmc, niter, nburnin))
```
```{r model6_obe_f7_densityplot, fig.height = 2, fig.align = "center",eval=TRUE}
                samplesPlot(samples, 'obe', densityplot = FALSE)
                samples_temp <- samples[1000:(niter-nburnin),]
                table<-samplesSummary(samples_temp)
                round(table,3) %>%
                kbl(caption = "Estimates") %>%
                kable_classic(full_width = F, html_font = "Cambria")  
```

## Support: OSE in presence of AE (model 6)

```{r model5_6_OSE_init, echo=FALSE}  
    N       <- nrow(df)
    sexPref <- df$sexPref
    os      <- df$os
    Nsib    <- df$Nsib
    yob     <- df$yob_s 
```

### OSE function 7

```{r model6_ose f7,echo=TRUE}
 niter <- 500000
 nburnin <- 10
 code <- nimbleCode({
           a1 ~ dunif(0,1)
      phi_ea ~ dunif(0, 1)
         phi ~ dunif(0, 1)
          p0 ~ dunif(0, 1)
       alpha ~ dunif(0 , 2)
        beta ~ dunif(0, 10)
        c[1] ~ dunif(-2, 10)
        c[2] ~ dunif(-10, 10)
        obe ~ dbern(0.5)
  beta_hyper ~ dunif(0, 10)
  for(i in 1:N) {
       lambda[i] <- exp(c[1] + c[2]*yob[i])*(1+beta*x[i])*(1+beta_hyper*x_hy[i])
       p1[i]<-min((1+alpha*x[i])*p0,0.999999)
         p[i]<- min(p1[i] + step(os[i]-1)*a1*obe,0.99999) 
         Nsib[i] ~ dpois(lambda[i])       
      sexPref[i] ~ dbern(p[i])
            x[i] ~ dbern(phi_ea)
         x_hy[i] ~ dbern(phi)
  }
})
constants = list(N=N,yob=yob,os=os)
data = list(Nsib=Nsib,sexPref=sexPref)
inits = list(a1=0.1, c=c(1,0.1), x = rbinom(N, size=1, prob=0.5),
              x_hy = rbinom(N, size=1, prob=0.5), beta_hyper = 1,obe=0,
              phi_ea = .5, phi=0.5,alpha = .5, beta = .5,p0=0.3)
```
```{r model6_ose_f7_suite, eval=TRUE}
        
        Rmodel <- nimbleModel(code, constants, data, inits)
        Rmodel$calculate()
        
        conf <- configureMCMC(Rmodel, enableWAIC = TRUE)
        conf$printSamplers(byType = TRUE)
        conf$removeSampler('obe')
        conf$printSamplers(byType = TRUE)
        conf$addSampler(target = 'obe', type = 'RJ_indicator1',   
                                control =
                                  list(coefNode1 = 'a1', propScale1 = 0.3))  
        conf$printSamplers(byType = TRUE)
        nodes <- c('a1')  
        print(nodes)
        
        for(node in nodes) {
          samplerConf <- conf$getSamplers(node)[[1]]
          conf$removeSamplers(node)
          conf$addSampler(target = node, type = 'RJ_toggle2',
                          control =
                            list(indicatorNode = 'obe',
                                 samplerConf = samplerConf))
        }
        conf$printSamplers(byType = TRUE)
        conf$printMonitors()
        Rmcmc <- buildMCMC(conf)
        Cmodel <- compileNimble(Rmodel)
        Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
        set.seed(round(runif(1, min=0, max=1000), digits=0))
        system.time(samples <- runMCMC(Cmcmc, niter, nburnin))
```
```{r model6_ose_f7_densityplot, fig.height = 2, fig.align = "center",eval=TRUE}
                samplesPlot(samples, 'obe', densityplot = FALSE)
                samples_temp <- samples[1000:(niter-nburnin),]
                table<-samplesSummary(samples_temp)
                round(table,3) %>%
                kbl(caption = "Estimates") %>%
                kable_classic(full_width = F, html_font = "Cambria")   #library(kableExtra)
```

## Support:  AE in presence of an OSE (model 5)

### OSE function 7

```{r model5_ose_f7, echo = TRUE}
niter <- 500000
nburnin <- 10
code <- nimbleCode({
           a1 ~ dunif(0,1)
      phi_ea ~ dunif(0, 1)
         phi ~ dunif(0, 1)
          p0 ~ dunif(0, 1)
       alpha ~ dunif(0 , 2)
        beta ~ dunif(0, 10)
        c[1] ~ dunif(-2, 10)
        c[2] ~ dunif(-10, 10)
        ant ~ dbern(0.5)
  beta_hyper ~ dunif(0, 10)
  for(i in 1:N) {
       lambda[i] <- exp(c[1] + c[2]*yob[i])*(1+beta*x[i]*ant)*(1+beta_hyper*x_hy[i])
       p1[i]<-min((1+alpha*x[i]*ant)*p0,0.999999)
       p[i]<- min(p1[i] + step(os[i]-1)*a1,0.99999) 
         Nsib[i] ~ dpois(lambda[i])       
      sexPref[i] ~ dbern(p[i])
            x[i] ~ dbern(phi_ea)
         x_hy[i] ~ dbern(phi)
  }
})
constants = list(N=N,yob=yob,os=os)
data = list(Nsib=Nsib,sexPref=sexPref)
inits = list(a1=0.1, c=c(1,0.1), x = rbinom(N, size=1, prob=0.5),
              x_hy = rbinom(N, size=1, prob=0.5), beta_hyper = 1,ant=0,
              phi_ea = .5, phi=0.5,alpha = .5, beta = .5,p0=0.3)

```
```{r model5_ose_f7_suite, eval=TRUE}
        
        Rmodel <- nimbleModel(code, constants, data, inits)
        Rmodel$calculate()
        conf <- configureMCMC(Rmodel, enableWAIC = TRUE)
        conf$printSamplers(byType = TRUE)
        conf$removeSampler('ant')
        conf$printSamplers(byType = TRUE)
        conf$addSampler(target = 'ant', type = 'RJ_indicator2',
                        control =
                          list(coefNode1 = 'alpha', propScale1 = 0.1,    
                               coefNode2 = 'beta',  propScale2 = 0.5)) 
        conf$printSamplers(byType = TRUE)
        nodes <- c('alpha', 'beta')  
        print(nodes)
        for(node in nodes) {
          samplerConf <- conf$getSamplers(node)[[1]]
          conf$removeSamplers(node)
          conf$addSampler(target = node, type = 'RJ_toggle2',
                          control =
                            list(indicatorNode = 'ant',
                                 samplerConf = samplerConf))
        }
        conf$printSamplers(byType = TRUE)
        conf$printMonitors()
        Rmcmc <- buildMCMC(conf)
        Cmodel <- compileNimble(Rmodel)
        Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
        set.seed(round(runif(1, min=0, max=1000), digits=0))
        system.time(samples <- runMCMC(Cmcmc, niter, nburnin))
```
```{r model5_ose_f7_densityplot, fig.height = 2, fig.align = "center",eval=TRUE}
                samplesPlot(samples, 'ant', densityplot = FALSE)
                samples_temp <- samples[1000:(niter-nburnin),]
                table<-samplesSummary(samples_temp)
                round(table,3) %>%
                kbl(caption = "Estimates") %>%
                kable_classic(full_width = F, html_font = "Cambria")  
```
