monogenea

# Try Eggs_laid ~ dpois
froReduced <- slice(fro, which(!is.na(Eggs_laid))) %>% 
      as.data.frame()

# Re-do the variable scaling, otherwise the sampling throws an error
froReduced %<>% mutate(female_id = as.integer(factor(Female_ID_coded)),
                       year_id = as.integer(factor(Year)),
                       group_id = as.integer(factor(Group_ID_coded)),
                       Min_age_Z = scale(Min_age),
                       Group_size_Z = scale(Group_size),

monogenea

eggsLMod <- map2stan(alist(
      Eggs_laid ~ dpois(lambda),
      log(lambda) <- a + a_fem[female_id] + a_year[year_id] + a_group[group_id] +
            Parasite*bP + Min_age_Z*bA + Group_size_Z*bGS + Mean_eggsize_Z*bES +
            Parasite*Min_age_Z*bPA,
      Group_size_Z ~ dnorm(0, 3),
      Mean_eggsize_Z ~ dnorm(0, 3),
      a_fem[female_id] ~ dnorm(0, sigma1),
      a_year[year_id] ~ dnorm(0, sigma2),
      a_group[group_id] ~ dnorm(0, sigma3),

monogenea

# Sample posterior
post <- extract.samples(eggsLMod)

# Run simulations w/ averages of all predictors, except parasite 0 / 1
lambdaNoP <- exp(post$a + 0*post$bP + 0*post$bA +
                       0*post$bGS + 0*post$bES + 0*0*post$bPA)
simFledgeNoPar <- rpois(n = length(lambdaNoP), lambda = lambdaNoP)

lambdaP <- exp(post$a + 1*post$bP + 0*post$bA +
                     0*post$bGS + 0*post$bES + 1*0*post$bPA)

monogenea

# Bonus! Sample counts from predictionsP, take 95% HDPI
hdpiPoisP <- apply(predictionsP, 2, HPDI, prob = .95)
meanPoisP <- colMeans(predictionsP)
plot(rangeA, meanPoisP, type = "l", ylim = c(0, 15),
     yaxp = c(0, 15, 5), xlab = "Min Age (standardized)",
     ylab = expression(paste(lambda, " / no. eggs laid")))
shade(hdpiPoisP, rangeA)
poisample <- sapply(1:100, function(k){
      rpois(nrow(predictionsP), predictionsP[,k])
})

monogenea

library(tensorflow)
use_condaenv("greta")
library(greta)
library(tidyverse)
library(bayesplot)
library(readxl)

# Read female reproductive output and discard records w/ NAs
fro <- read_xlsx("data.xlsx", sheet = allTabs[2])
fro <- fro[complete.cases(fro),]

monogenea

# Define model effects
sigmaML <- cauchy(0, 1, truncation = c(0, Inf), dim = 3)
a_fem <- normal(0, sigmaML[1], dim = max(female_id))
a_year <- normal(0, sigmaML[2], dim = max(year))
a_group <- normal(0, sigmaML[3], dim = max(group_id))

a <- normal(0, 5)
bA <- normal(0, 3)
bEL <- normal(0, 3)
bES <- normal(0, 3)

monogenea

# Simulation with average eggs laid, egg size and group size, w/ and w/o parasitism
seqX <- seq(-3, 3, length.out = 100)
probsNoPar <- sapply(seqX, function(x){
      scenario <- ilogit(a + x*bA)
      probs <- calculate(scenario, draws)
      return(unlist(probs))
})
probsPar <- sapply(seqX, function(x){
      scenario <- ilogit(a + x*bA + bP + x*bPA)
      probs <- calculate(scenario, draws)

monogenea

# Mon Oct 29 13:17:55 2018 ------------------------------
## DREAM olfaction prediction challenge
library(caret)
library(rsample)
library(tidyverse)
library(recipes)
library(magrittr)
library(doMC)

# Create directory and download files

monogenea

# Determine intersection of compounds in features and responses
commonMols <- intersect(responses$CID,
                        molFeats$CID)
# Subset features and responses accordingly
responses %<>% filter(CID %in% commonMols)
molFeats %<>% filter(CID %in% commonMols)

# Compute median pleasantness across the population
medianPlsnt <- responses %>% 
      group_by(CID) %>% 

monogenea

# Filter nzv
X <- X[,-nearZeroVar(X, freqCut = 4)] # == 80/20