In famuvie/2016_RisingAshes: Companion to the Paper "Rising Out of the Ashes"

## Required packages
library(RisingAshes)
library(tidyr)    # data manipulation
library(plyr)     # data manipulation
library(dplyr)    # data manipulation
library(ggplot2)  # plotting
library(grid)     # plotting
library(gridExtra)# plotting
library(viridis)  # color palette
library(knitr)    # kable()
library(INLA)     # Bayesian Inference

## knitr options
opts_chunk$set(echo       = FALSE,
               message    = FALSE,
               warning    = FALSE,
               comment    = NA,
               fig.dpi    = 96,
               fig.width  = 5,
               fig.height = 3,
               fig.caption= FALSE,
               cache      = FALSE)

## ggplot2 options
theme_set(theme_bw())

## Applies an square-root transform to CD (creating variables CDT)
## Buils object mdat in long format (with variables: Year; Variable; Value)
## Open with:
# file.show(file=system.file('vignettes/cleanup.R',
#                            package = 'RisingAshes'))
source('cleanup.R')

# Object for storing models and solutions
ml <- list()

## Base dataset for analyzing the CL variable
ml$dat <- 
  mdat %>% 
  filter(
    variable %in% c('CL', 'BC')
    ,Year %in% c('2012', '2013')
  ) %>%
  spread(variable, value) %>% 
  transform(BLC = factor(BLC),
            PROVxBLC = factor(PROV:factor(BLC)),
            FAMxBLC = factor(FAM:factor(BLC))) %>%
  droplevels()

These data are characterized by a large proportion of zeros. This calls for a mixture model with a probability $p_{ij}$ of observing the disease (i.e. $CL > 0$), and a continuous distribution for the observations greater than 0.

ggplot(ml$dat, aes(x = CL)) + geom_histogram(bins = 30)

# # Only the non-null data
# ggplot(data.frame(CL = ml$dat$CL[ml$dat$CL > 0]), aes(x = CL)) + geom_histogram()
# # Its log-transform (not very Gaussian...)
# ggplot(data.frame(CL = log(ml$dat$CL[ml$dat$CL > 0])), aes(x = CL)) + geom_histogram()
# tf <- function(x) log(x) - log(1 - x + 0.02)
# ggplot(data.frame(CL = tf(ml$dat$CL[ml$dat$CL > 0])), aes(x = CL)) + geom_histogram()
# # Its logit-transform (not very Gaussian...)
# ggplot(data.frame(CL = logit(ml$dat$CL[ml$dat$CL > 0])), aes(x = CL)) + geom_histogram()

ggplot(ml$dat, aes(X, Y)) +
  facet_wrap(~ Year) +
  geom_tile(aes(fill = CL, color = CL)) +
  coord_fixed() +
  scale_fill_viridis() +
  scale_color_viridis()

## Creates objects:
##   gen_str
##     $Ainv
##     $recoding
##     $contraint.matrix
##   prec.A
## Open with:
# file.show(file=system.file('vignettes/genetic_structure.R',
#                            package = 'RisingAshes'))
source('genetic_structure.R')

## Generates an object
## sp_str
##   $loc
##   $mesh
##   $spde
##   $priors
sp_str <- Devecey_spde()

## mesh space (mesh$n * nYears = 942 * 5 = 4710 rows)
st.index <- inla.spde.make.index(
  name = 'theta',
  n.spde = sp_str$spde$n.spde,
  n.group = nlevels(ml$dat$Year)
)

## Maps from the observation to the prediction space
st.Amat <- inla.spde.make.A(
  sp_str$mesh, loc = sp_str$loc, 
  index = rep(1:nrow(sp_str$loc), nlevels(ml$dat$Year)),
  group = as.numeric(ml$dat$Year)
)

Models for the binary component

Role of the Basal Circumference

Here we compare several models to take into account the Basal Circumference. It is quite clear that the best model is the one with categories interacting with the Year.

There is a strong jump in the behaviour of resistance for trees under 30 cm of circumference, with respect to those above. However, there is also some variation in the other categories as well.

## Models for the binary component
## with Basal Circumference as explanatory variable

ml$bin_BC$dat <- 
  ml$dat[, c('Seq', 'FAM', 'PROV', 'BF98', 'Year', 'BC')] %>% 
  transform(
    Year2013 = as.numeric(Year == 2013),
    mu = 1,
    BC.cat = cut(BC, c(0, 30, 45, 60, 75, 90, 150), ordered = TRUE),
    BC.sc = scale(BC),
    recode = gen_str$recoding,
    FAM = as.numeric(FAM),
    CL0 = as.numeric(ml$dat$CL == 0)
  )

## Binomial stack
ml$bin_BC$stk <-
  inla.stack(
    data    = list(y = ml$bin_BC$dat$CL0),
    A       = list(st.Amat, 1),
    effects = list(
      suffix(st.index, 'bin'),
      suffix(ml$bin_BC$dat, 'bin')
    ),
    tag     = 'bin'
  )

## Base Formula: without BC
ml$bin_BC$fml <- 
  y ~ 0 + Year.bin + 
  f(theta.bin, model=sp_str$spde,
    group = theta.group.bin,  # variable in the mesh space
    control.group = list(model = 'exchangeable')) +
  f(recode.bin, model = "generic0", hyper = list(theta = prec.A),
    constr = TRUE, Cmatrix = Cmat)

## Fit models for comparison
## Open with:
# file.show(file=system.file('vignettes/CL_bin_vs_BC_models.R',
#                            package = 'RisingAshes'))
source('CL_bin_vs_BC_models.R')

## Retrieve results from server
ml$bin_BC$res <- c(lapply(ml$bin_BC$res, inla.qget))

## Load into ml$bin_BC the cached version of the two previous chunks
ml$bin_BC <- RisingAshes:::BCmodels_CLbin

lapply(ml$bin_BC$res, summary)

binBC_mcdat <- 
  data.frame(
    DIC = sapply(ml$bin_BC$res,
                 function(x) x$dic$dic),
    p_D = sapply(ml$bin_BC$res,
                 function(x) x$dic$p.eff),
    WAIC = sapply(ml$bin_BC$res,
                  function(x) x$waic$waic),
    p_W = sapply(ml$bin_BC$res,
                 function(x) x$waic$p.eff),
    marginal_loglik = sapply(ml$bin_BC$res,
                             function(x) x$mlik[1, 1]),
    LPML = sapply(ml$bin_BC$res,
                  function(x) sum(log(x$cpo$cpo), na.rm = TRUE))
  )

kable(
  binBC_mcdat,
  digits = 0
)

ggplot(binBC_mcdat[-1, ], aes(DIC, marginal_loglik)) +
  geom_point() + 
  geom_text(aes(label = rownames(binBC_mcdat)[-1]), hjust = -.1)

BCdat <- local({
  ## Values of Basal Circumference to evaluate predictions
  x <- seq(min(ml$bin_BC$dat$BC, na.rm = TRUE),
           max(ml$bin_BC$dat$BC, na.rm = TRUE),
           length = 51)


  ## Predictions of the categorical model in the values of x
  x.cat <-  cut(x, c(0, 30, 45, 60, 75, 90, 150), ordered = TRUE)

  cat.values <- 
    with(
      ml$bin_BC$res$cat,
      c(summary.fixed[1:6, 'mean'],
        summary.fixed[7, 'mean'] + summary.fixed[1:6, 'mean'])
    )

  cat <- c(cat.values[x.cat], cat.values[7:12][x.cat])

  catxyear.values <- 
    with(
      ml$bin_BC$res$catxyear,
      as.vector(t(matrix(summary.fixed[, 'mean'], nrow = 2)))
    )

  catxyear <- c(catxyear.values[x.cat], catxyear.values[7:12][x.cat])

  cat2.values <-
    with(ml$bin_BC$res$cat2,
         as.vector(t(matrix(summary.fixed[, 'mean'], nrow = 2)))
    )

  cat2 <- c(cat2.values[cut(x, c(0, 30, 150))],
            cat2.values[3:4][cut(x, c(0, 30, 150))])

  ## Predictions of the linear and quadratic models in the values of x
  x.lin <- 
    with(
      attributes(scale(ml$bin_BC$dat$BC)),
      scale(x,
            center = `scaled:center`,
            scale = `scaled:scale`)
    )

  lin <- 
    with(
      ml$bin_BC$res$lin,
      c(summary.fixed['Year.bin2012', 'mean'] + 
          summary.fixed['BC.sc.bin', 'mean'] * x.lin,
        summary.fixed['Year.bin2013', 'mean'] + 
          summary.fixed['BC.sc.bin', 'mean'] * x.lin)
    )

  quad <- 
    with(
      ml$bin_BC$res$quad,
      c(summary.fixed['Year.bin2012', 'mean'] + 
          summary.fixed['BC.sc.bin', 'mean'] * x.lin +
          summary.fixed['I(BC.sc.bin^2)', 'mean'] * x.lin^2,
        summary.fixed['Year.bin2013', 'mean'] + 
          summary.fixed['BC.sc.bin', 'mean'] * x.lin +
          summary.fixed['I(BC.sc.bin^2)', 'mean'] * x.lin^2)
    )

  np  <-
    with(
      ml$bin_BC$res$np,
      c(summary.fixed['Year.bin2012', 'mean'] +
          summary.random$`cut(BC.sc.bin, 51)`[, 'mean'],
        summary.fixed['Year.bin2013', 'mean'] +
          summary.random$`cut(BC.sc.bin, 51)`[, 'mean'])
    )

  data.frame(
    Year = factor(rep(2012:2013, each = 51)),
    BC   = rep(x, 2),
    cat,
    cat2,
    catxyear,
    lin,
    quad,
    np) %>% 
    gather(model, y, cat:np)
})

ggplot(BCdat, aes(BC, y, colour = model)) +
  geom_line() + 
  facet_wrap(~ Year)

## Back-transform the posterior means, into the prob. scale
## Note this is not the posterior mean in that scale, but approx.

invlogit <- function(x) exp(x)/(1+exp(x))


tdat <- filter(mdat,
               variable %in% c('CD', 'CL', 'BC')
               ,Year %in% c('2012', '2013')
) %>%
  spread(variable, value) %>% 
  transform(CD = as.ordered(CD),
            BC = cut(BC, c(0, 30, 45, 60, 75, 90, 150), ordered = TRUE),
            CLbin = factor(CL > 0, labels = c('0', '>0'))) %>%
  filter(!is.na(CL), !is.na(CD)) %>% 
  droplevels() %>% 
  group_by(Year, BC) %>% 
  summarise(NCL0 = sum(CL == 0),
            TCL = n()) %>% 
  mutate(y = NCL0/TCL,
         PT.se = sqrt(y*(1-y)/TCL),
         BC = (15 * (.5 + 1:6))[BC])

ggplot(mutate(BCdat, y = invlogit(y)), aes(BC, y)) +
  geom_line(aes(colour = model)) + 
  geom_point(dat = tdat, stat = 'identity') +
  facet_wrap(~ Year) +
  ylab("Probability")

In conclusion, the BC variable seems to explain part of the binary component of CL. The main conclusion being that smaller trees have higher chance of escaping or resisting the disease.

However, BC is in turn explained in part by local environmental conditions. I.e. trees grow less on dryer areas, which are coincidentally bad conditions for the spread of the fungi. So the question arises to whether bad environmental condition causes less growth which in turn increases the chance of resisting the disease, or in reality bad environmental conditions causes both less growth and and increase in escaping the disease.

In fact, we see that while the models including the BC variable have lower DICs, the model without the variable have higher likelihood. This can be a result of a lack of separability between the spatial effect and the BC covariate.

Sice it is hard to find a biological relationship between BC and CL escaping/resisting, we stick to the theory that the local environmental conditions are the underlying cause of the apparent relationship. In consequence, we do not include this variable in the model.

Role of the Provenance

Here we test the relevance of the Provenance in explaining the Binary component.

While the DIC suggests using the main effect of the provenance of the covariate, the marginal likelihood of the model without provenance is somewhat larger. Furthermore, the sizes of the effects are very small. Thus, we consider that if there is some effect of the provenance then it must be very weak, and we choose not to use it in the model.

ml$bin_PROV <- ml$bin_BC[c('dat', 'stk', 'fml')]

## Run models for comparison
source('CL_bin_vs_PROV_models.R')

## Retrieve results from server
ml$bin_PROV$res <- c(lapply(ml$bin_PROV$res, inla.qget))

## Summary of spatial results
ml$bin_PROV$spderes <- 
  lapply(ml$bin_PROV$res[-1], inla.spde2.result, 'theta.bin', sp_str$spde)

## Load into ml$bin_PROV the cached version of the previous chunks
ml$bin_PROV <- RisingAshes:::PROVmodels_CLbin

lapply(ml$bin_PROV$res, summary)

binPROV_mcdat <- 
    data.frame(
    DIC = sapply(ml$bin_PROV$res,
                 function(x) x$dic$dic),
    p_D = sapply(ml$bin_PROV$res,
                 function(x) x$dic$p.eff),
    WAIC = sapply(ml$bin_PROV$res,
                  function(x) x$waic$waic),
    p_W = sapply(ml$bin_PROV$res,
                 function(x) x$waic$p.eff),
    marginal_loglik = sapply(ml$bin_PROV$res,
                      function(x) x$mlik[1, 1]),
    LPML = sapply(ml$bin_PROV$res,
                      function(x) sum(log(x$cpo$cpo), na.rm = TRUE))
  )

kable(
  binPROV_mcdat,
  digits = 0
)

ggplot(binPROV_mcdat[-1, ], aes(DIC, marginal_loglik)) +
  geom_point() + 
  geom_text(aes(label = rownames(binPROV_mcdat)[-1]), hjust = -.1)

PROVdat <- local({

  ## Predictions of the categorical model

  cat <- with(
    ml$bin_PROV$res$cat,
    rbind(summary.fixed[1:3, c('mean', '0.025quant', '0.975quant')],
          summary.fixed[4, 'mean'] +
            summary.fixed[1:3, c('mean', '0.025quant', '0.975quant')])
  )


  catxyear <- with(
    ml$bin_PROV$res$catxyear,
        summary.fixed[t(matrix(1:6, 2)),
                      c('mean', '0.025quant', '0.975quant')]
  )


  ## mean of the sum of ST and genetic random effects
  ## for 2012 and 2013, in the cat model.
  center.ranef <-
    sapply(
    1:2,
    function(x) 
      mean(ml$bin_PROV$spderes$cat$summary.values$mean[st.index$theta.group == x],
           na.rm = TRUE)) +
    mean(ml$bin_PROV$res$cat$summary.random$recode.bin$mean[gen_str$recoding], 
         na.rm = TRUE)



  data.frame(
    Year = factor(rep(2012:2013, each = 3)),
    centering = rep(center.ranef, each = 3),
    PROV   = factor(rep(levels(ml$dat$PROV), 2)),
    cat = cat,
    cy = catxyear ) %>% 
    gather(model, y, cat.mean:cy.0.975quant) %>% 
    separate(model, c('model', 'q'), extra = 'merge') %>% 
    mutate(y = y + centering) %>% 
    spread(q, y)
})

ggplot(PROVdat, aes(PROV, mean, colour = model)) +
  geom_pointrange(aes(ymin = `0.025quant`, ymax = `0.975quant`),
                  position = position_dodge(width = 0.5)) +
  facet_wrap(~ Year)

## Back-transform the posterior means, into the prob. scale
## Note this is not the posterior mean in that scale, but approx.

invlogit <- function(x) exp(x)/(1+exp(x))


tdat <- filter(mdat,
               variable %in% c('CD', 'CL')
               ,Year %in% c('2012', '2013')
) %>%
  spread(variable, value) %>% 
  transform(CLbin = factor(CL > 0, labels = c('0', '>0'))) %>%
  filter(!is.na(CL), !is.na(PROV)) %>% 
  droplevels() %>% 
  group_by(Year, PROV) %>% 
  summarise(NCL0 = sum(CL == 0),
            TCL = n()) %>% 
  mutate(y = NCL0/TCL,
         PT.se = sqrt(y*(1-y)/TCL))

ggplot(mutate(PROVdat,
              y = invlogit(mean)), aes(PROV, y)) +
  geom_bar(aes(fill = model), stat = 'identity', position = 'dodge') + 
  geom_point(dat = tdat, stat = 'identity') +
  facet_wrap(~ Year) +
  ylab("Probability")

Role of the Bud-flush

Here we test the relevance of the Bud-flush in explaining the Binary component.

ml$bin_BF <- ml$bin_BC[c('dat', 'stk', 'fml')]

## Run models for comparison
source('CL_bin_vs_BF_models.R')

ml$bin_BF$res <- c(lapply(ml$bin_BF$res, inla.qget))

## Summary of spatial results
ml$bin_BF$spderes <- lapply(ml$bin_BF$res[-1], inla.spde2.result, 'theta.bin', sp_str$spde)

## Load into ml$bin_BF the cached version of the previous chunks
ml$bin_BF <- RisingAshes:::BFmodels_CLbin

lapply(ml$bin_BF$res, summary)

BF_mcdat <- 
    data.frame(
    DIC = sapply(ml$bin_BF$res,
                 function(x) x$dic$dic),
    WAIC = sapply(ml$bin_BF$res,
                  function(x) x$waic$waic),
    marginal_loglik = sapply(ml$bin_BF$res,
                      function(x) x$mlik[1, 1]),
    LPML = sapply(ml$bin_BF$res,
                      function(x) sum(log(x$cpo$cpo), na.rm = TRUE))
  )

kable(
  BF_mcdat,
  digits = 0
)

ggplot(BF_mcdat[-1, ], aes(DIC, marginal_loglik)) +
  geom_point() + 
  geom_text(aes(label = rownames(BF_mcdat)[-1]), hjust = -.1)

BFdat <- local({

  nlevels <- 5
  ## Predictions of the categorical model

  cat <- with(
    ml$bin_BF$res$cat,
    rbind(summary.fixed[1:nlevels, c('mean', '0.025quant', '0.975quant')],
          summary.fixed[nlevels+1, 'mean'] +
            summary.fixed[1:nlevels, c('mean', '0.025quant', '0.975quant')])
  )


  catxyear <- with(
    ml$bin_BF$res$catxyear,
        summary.fixed[t(matrix(1:(2*nlevels), 2)),
                      c('mean', '0.025quant', '0.975quant')]
  )


  center.ranef <-
    sapply(
      1:2,
      function(x) 
        mean(ml$bin_BF$spderes$cat$summary.values$mean[st.index$theta.group == x],
             na.rm = TRUE)) +
    mean(ml$bin_BF$res$cat$summary.random$recode.bin$mean[gen_str$recoding],
         na.rm = TRUE)



  data.frame(
    Year = factor(rep(2012:2013, each = nlevels)),
    centering = rep(center.ranef, each = nlevels),
    BF   = factor(rep(levels(ml$dat$BF98), 2)),
    cat = cat,
    cy = catxyear ) %>% 
    gather(model, y, cat.mean:cy.0.975quant) %>% 
    separate(model, c('model', 'q'), extra = 'merge') %>% 
    mutate(y = y + centering) %>% 
    spread(q, y)
})

ggplot(BFdat, aes(BF, mean, colour = model)) +
  geom_pointrange(aes(ymin = `0.025quant`, ymax = `0.975quant`),
                  position = position_dodge(width = 0.5)) +
  facet_wrap(~ Year)

## Back-transform the posterior means, into the prob. scale
## Note this is not the posterior mean in that scale, but approx.

invlogit <- function(x) exp(x)/(1+exp(x))


tdat <- filter(mdat,
               variable %in% c('CD', 'CL')
               ,Year %in% c('2012', '2013')
) %>%
  spread(variable, value) %>% 
  transform(CLbin = factor(CL > 0, labels = c('0', '>0'))) %>%
  filter(!is.na(CL), !is.na(BF98)) %>% 
  droplevels() %>% 
  group_by(Year, BF98) %>% 
  summarise(NCL0 = sum(CL == 0),
            TCL = n()) %>% 
  mutate(y = NCL0/TCL,
         PT.se = sqrt(y*(1-y)/TCL),
         BF = BF98)

ggplot(mutate(BFdat,
              y = invlogit(mean)), aes(BF, y)) +
  geom_bar(aes(fill = model), stat = 'identity', position = 'dodge') + 
  geom_point(dat = tdat, stat = 'identity') +
  facet_wrap(~ Year) +
  ylab("Probability")

As a conclusion, while there seems to be a trend in the relationship between the binary component and the BF, it is not strong enough to justify contemplating a variable that is more naturally linked to the crown than to the base. Therefore, we choose not to include this variable in the model.

Model selection binary component

We have analysed the effects of the Basal Circumference (BC), the Provenance (PROV) and the Bud-flushing (BF) separately. In all three cases the estimated effects are small and the marginal likelihood is highest for the model without the variable included in any form. On the other hand, the other scores are optimal on one of the other models.

In the case of BC we hypothesize a confusion with the spatial effect. The PROV is clearly not relevant. Finally, the BF seems to be somewhat related, but it is hard to justify biologically.

Based on this measures and on the effects sizes we make a preliminary variable selection as follows:

include the BC variable in its categorical form and its interaction with the year
do not include the PROV variable
do not consider the BF variable

Now we compare these three selected submodels (only the one with BC do not coincide with the null), with the null model (not including any of those variables) and with the full model (including both BC and BF).

## Models for the binary component

ml$bin_msel$dat <- 
  ml$dat[, c('Seq', 'FAM', 'PROV', 'BF98', 'Year', 'BC')] %>% 
  transform(
    Year2013 = as.numeric(Year == 2013),
    mu = 1,
    BC = cut(BC, c(0, 30, 45, 60, 75, 90, 150), ordered = TRUE),
    BF_2 = as.numeric(BF98 == 2),
    BF_3 = as.numeric(BF98 == 3),
    BF_4 = as.numeric(BF98 == 4),
    BF_5 = as.numeric(BF98 == 5),
    recode = gen_str$recoding,
    FAM = as.numeric(FAM),
    CL0 = as.numeric(ml$dat$CL == 0))

## Binomial stack
ml$bin_msel$stk <-
  inla.stack(data    = list(y = ml$bin_msel$dat$CL0),
             A       = list(st.Amat, 1),
             effects = list(suffix(st.index, 'bin'),
                            suffix(ml$bin_msel$dat, 'bin')),
             tag     = 'bin')

## Base Formula: without BC
ml$bin_msel$fml <- y ~ 0 + Year.bin + 
               f(theta.bin, model=sp_str$spde,
                 group = theta.group.bin,  # variable in the mesh space
                 control.group = list(model = 'exchangeable')) +
               f(recode.bin, model = "generic0", hyper = list(theta = prec.A),
                 constr = TRUE, Cmatrix = Cmat)

## Run the full model
source('CL_bin_models.R')

ml$bin_msel$res <- c(lapply(ml$bin_msel$res, inla.qget))

## Load into ml$bin_msel$res$full the cached version of the two previous chunks
ml$bin_msel$res$full <- RisingAshes:::full_CLbin

## These were already fitted
ml$bin_msel$res <- c(list(
  null = ml$bin_BC$res$noBC,
  BC   = ml$bin_BC$res$catxyear,
  BF   = ml$bin_BC$res$cat
),
ml$bin_msel$res
)

## It would be good to fit also all the combinations with two of those

lapply(ml$bin_msel$res, summary)

BF_mcdat <- 
    data.frame(
    DIC = sapply(ml$bin_msel$res,
                 function(x) x$dic$dic),
    WAIC = sapply(ml$bin_msel$res,
                  function(x) x$waic$waic),
    marginal_loglik = sapply(ml$bin_msel$res,
                      function(x) x$mlik[1, 1]),
    LPML = sapply(ml$bin_msel$res,
                      function(x) sum(log(x$cpo$cpo), na.rm = TRUE))
  )

kable(
  BF_mcdat,
  digits = 0
)

ggplot(BF_mcdat, aes(DIC, marginal_loglik)) +
  geom_point() + 
  geom_text(aes(label = rownames(BF_mcdat)), hjust = -.1)

We find an antagonistic relationship between two selection criteria. Based on previous observations, we suspect that this is an effect of a lack of identifiability in the model. Therefore, we choose not to include these explanatory variables in the binary component of the model.

Models for the continuous component

We perform a preliminary likelihood evaluation, to assess whether the data is better fitted using a Beta or a Gamma likelihood.

Both optios have their limitations. The Beta distribution does not admit values of 1 (that are arbitrarily switched to 0.98 here). On the other hand, the Gamma has a semi-infinite support (while we know that there is a maximum at 1 which can eventually accumulate probability mass).

ml$cont$dat <- 
  ml$dat[, c('Seq', 'FAM', 'PROV', 'BF98', 'Year', 'BC')] %>% 
  transform(
    Year2013 = as.numeric(Year == 2013),
    BC.cat = cut(BC, c(0, 30, 45, 60, 75, 90, 150), ordered = TRUE),
    BC.sc = scale(BC),
    BF_2 = as.numeric(BF98 == 2),
    BF_3 = as.numeric(BF98 == 3),
    BF_4 = as.numeric(BF98 == 4),
    BF_5 = as.numeric(BF98 == 5),
    recode = gen_str$recoding,
    CL0 = as.numeric(ml$dat$CL == 0))

ml$cont$dat$CL_gt0          <- ml$dat$CL
ml$cont$dat$CL_gt0[ml$dat$CL == 0] <- NA

ml$cont$dat$yt <- ml$cont$dat$CL_gt0
ml$cont$dat$yt[ml$dat$CL == 1] <- .98   # forbiden value for a beta

## Use the same covariates as before

fit_lik <- function(lik) {
  fml <- yt ~ 0 + FAM + PROV + Year:BC
  inla(formula = fml,
       data = filter(ml$cont$dat,
                     !is.na(BF98),
                     !is.na(FAM)),
       family = lik,
       control.predictor = list(
         compute = TRUE, link = 1),
       control.fixed = list(
         expand.factor.strategy = 'inla'), # account for NAs in factor
       control.compute = list(
         dic = TRUE,
         waic = TRUE,
         cpo = TRUE,
         po = TRUE,
         config = TRUE)
       # , inla.call = 'remote'
  )
}

liklst <- c('beta', 'gamma')
cont_res <- lapply(liklst, fit_lik)
names(cont_res) <- liklst

cont_comp <- 
    data.frame(
    DIC = sapply(cont_res,
                 function(x) x$dic$dic),
    p_D = sapply(cont_res,
                 function(x) x$dic$p.eff),
    WAIC = sapply(cont_res,
                  function(x) x$waic$waic),
    p_W = sapply(cont_res,
                 function(x) x$waic$p.eff),
    marginal_loglik = sapply(cont_res,
                      function(x) x$mlik[1, 1]),
    LPML = sapply(cont_res,
                      function(x) sum(log(x$cpo$cpo), na.rm = TRUE))
  )

kable(
  cont_comp,
  digits = 0
)

The results clearly favour the Gamma distribution.

In the following models we include an overdispersion term in the latent space. This is simply unstructured noise, with mean zero and variance $\sigma^2_{\text{od}}$. On exchange, we fix a very high precision for the likelihood distribution.

In this way, we include the distribution-specific variance in the observations scale by an additional random component in the latent scale. This equivalent formulation helps in the specification of a prior distribution for its variance, and facilitates the computation of the heritability.

On the downside, using this trick we loose the ability to compare the models using DIC or WAIC. They are all the same as a consequence of fixing the precision of the likelihood.

Update: step-back from this strategy. I can save it for the final model. For model comparison is more convenient to have access to all measures of model fit. Even if the prior on the precision of the likelihood is approximate, it is the same for all models, so that the comparison is fair. On the other hand, it is not that difficult to assess. The scale should be also around the unit. Therefore, we can safely use the same prior that we were using for the precision of the latent overdispersion parameter.

## Models for the continuous component
## with Basal Circumference as explanatory variable


## continuous stack

ml$cont$stk <-
  inla.stack(
    data    = list(y = ml$cont$dat$yt),
    A       = list(st.Amat, 1),
    effects = list(suffix(st.index, 'cont'),
                   suffix(ml$cont$dat, 'cont')),
    tag     = 'cont')


## Prior on the precision of the overdispersion term
## We use a Gamma with shape 0.5 and rate 0.2, which gives a
## prior expectation (in the log scale) of 2.5 and a reasonably ample
## variance of 12.5.
## Furthermore, we checked that sensitivity to this prior specification is very
## limited.
prec.od <- list(fixed = FALSE,
                initial = 1,
                param = c(.5, .2))


## Base Formula: without any covariate but Year
ml$cont$fml <- 
  y ~ 0 + Year.cont + 
  f(theta.cont, model=sp_str$spde,
    group = theta.group.cont,  # variable in the mesh space
    control.group = list(model = 'exchangeable')) +
  f(recode.cont, model = "generic0", hyper = list(theta = prec.A),
    constr = TRUE, Cmatrix = Cmat) #+
  # f(Seq.cont, model = 'iid', hyper = list(theta = prec.od))

Role of the Basal Circumference

Here we compare several models to take into account the Basal Circumference. It is quite clear that the best model is the one with categories interacting with the Year.

There is a strong jump in the behaviour of resistance for trees under 30 cm of circumference, with respect to those above. However, there is also some variation in the other categories as well.

## Run models for comparison
source('CL_cont_vs_BC_models.R')

ml$cont_BC$res <- c(lapply(ml$cont_BC$res, inla.qget))

## Load into ml$cont_BC the cached version of the previous chunks
ml$cont_BC <- RisingAshes:::BCmodels_CLcont

lapply(ml$cont_BC$res, summary)

contBC_mcdat <- 
    data.frame(
    DIC = sapply(ml$cont_BC$res,
                 function(x) x$dic$dic),
    p_D = sapply(ml$cont_BC$res,
                 function(x) x$dic$p.eff),
    WAIC = sapply(ml$cont_BC$res,
                  function(x) x$waic$waic),
    p_W = sapply(ml$cont_BC$res,
                 function(x) x$waic$p.eff),
    marginal_loglik = sapply(ml$cont_BC$res,
                      function(x) x$mlik[1, 1]),
    LPML = sapply(ml$cont_BC$res,
                      function(x) sum(log(x$cpo$cpo), na.rm = TRUE))
  )

kable(
  contBC_mcdat,
  digits = 0
)

ggplot(contBC_mcdat[-1, ], aes(DIC, marginal_loglik)) +
  geom_point() + 
  geom_text(aes(label = rownames(contBC_mcdat)[-1]), hjust = -.1)

BCdat <- local({
  ## Values of Basal Circumference to evaluate predictions
  x <- seq(min(ml$cont$dat$BC, na.rm = TRUE),
           max(ml$cont$dat$BC, na.rm = TRUE), 
           length = 51)


  ## Predictions of the categorical model in the values of x
  x.cat <-  cut(x, c(0, 30, 45, 60, 75, 90, 150), ordered = TRUE)

  cat.values = with(ml$cont_BC$res$cat,
                 c(summary.fixed[1:6, 'mean'],
                   summary.fixed[7, 'mean'] + summary.fixed[1:6, 'mean'])
  )

  cat <- c(cat.values[x.cat], cat.values[7:12][x.cat])


  catxyear.values = with(ml$cont_BC$res$catxyear,
                    as.vector(t(matrix(summary.fixed[, 'mean'], nrow = 2)))
  )

  catxyear <- c(catxyear.values[x.cat], catxyear.values[7:12][x.cat])

  cat2.values = with(ml$cont_BC$res$cat2,
                         as.vector(t(matrix(summary.fixed[, 'mean'], nrow = 2)))
  )

  cat2 <- c(cat2.values[cut(x, c(0, 30, 150))], cat2.values[3:4][cut(x, c(0, 30, 150))])


  ## Predictions of the linear and quadratic models in the values of x
  x.lin <- with(
    attributes(scale(ml$cont$dat$BC)),
    scale(x,
          center = `scaled:center`,
          scale = `scaled:scale`)
  )

  lin = with(
    ml$cont_BC$res$lin,
    c(summary.fixed['Year.cont2012', 'mean'] + 
        summary.fixed['BC.sc.cont', 'mean'] * x.lin,
      summary.fixed['Year.cont2013', 'mean'] + 
        summary.fixed['BC.sc.cont', 'mean'] * x.lin)
  )

  quad = with(
    ml$cont_BC$res$quad,
    c(summary.fixed['Year.cont2012', 'mean'] + 
        summary.fixed['BC.sc.cont', 'mean'] * x.lin +
        summary.fixed['I(BC.sc.cont^2)', 'mean'] * x.lin^2,
      summary.fixed['Year.cont2013', 'mean'] + 
        summary.fixed['BC.sc.cont', 'mean'] * x.lin +
        summary.fixed['I(BC.sc.cont^2)', 'mean'] * x.lin^2)
  )


  np = with(
    ml$cont_BC$res$np,
    c(summary.fixed['Year.cont2012', 'mean'] +
        summary.random$`cut(BC.sc.cont, 51)`[, 'mean'],
      summary.fixed['Year.cont2013', 'mean'] +
        summary.random$`cut(BC.sc.cont, 51)`[, 'mean'])
  )


  data.frame(
    Year = factor(rep(2012:2013, each = 51)),
    BC   = rep(x, 2),
    cat,
    cat2,
    catxyear,
    lin,
    quad,
    np) %>% 
    gather(model, y, cat:np)
})

ggplot(BCdat, aes(BC, y, colour = model)) +
  geom_line() + 
  facet_wrap(~ Year)

## Back-transform the posterior means, into the prob. scale
## Note this is not the posterior mean in that scale, but approx.

tdat <- 
  filter(mdat,
         variable %in% c('CD', 'CL', 'BC')
         ,Year %in% c('2012', '2013')
  ) %>%
  spread(variable, value) %>% 
  transform(CD = as.ordered(CD),
            BC = cut(BC, c(0, 30, 45, 60, 75, 90, 150), ordered = TRUE),
            CLcont = factor(CL > 0, labels = c('0', '>0'))) %>%
  filter(!is.na(CL), !is.na(CD)) %>% 
  droplevels() %>% 
  group_by(Year, BC) %>% 
  summarise(y = mean(CL[CL > 0]),
            n = n()) %>% 
  mutate(BC = (15 * (.5 + 1:6))[BC])

ggplot(mutate(BCdat, y = exp(y)), aes(BC, y)) +
  geom_line(aes(colour = model)) + 
  geom_point(aes(size = n), dat = tdat, stat = 'identity') +
  facet_wrap(~ Year) +
  ylab("CL | CL > 0")

According to the DIC, the cat2 specification (BC in two categories interacting with the Year) is only slightly better than the model without BC, which has higher likelihood. All the other models using BC are far behind.

It looks that there is no real effect of BC, and that the seemingly better performance of cat2 is due to a very influential observation in 2012 with BC of 30 cm.

Role of the Provenance

Here we test the relevance of the Provenance in explaining the Binary component.

## Run models for comparison
source('CL_cont_vs_PROV_models.R')

ml$cont_PROV$res <- c(lapply(ml$cont_PROV$res, inla.qget))

## Summary of spatial results
ml$cont_PROV$spderes <- lapply(ml$cont_PROV$res[-1], inla.spde2.result, 'theta.cont', sp_str$spde)

## Load into ml$cont_PROV the cached version of the previous chunks
ml$cont_PROV <- RisingAshes:::PROVmodels_CLcont

lapply(ml$cont_PROV$res, summary)

contPROV_mcdat <- 
    data.frame(
    DIC = sapply(ml$cont_PROV$res,
                 function(x) x$dic$dic),
    p_D = sapply(ml$cont_PROV$res,
                 function(x) x$dic$p.eff),
    WAIC = sapply(ml$cont_PROV$res,
                  function(x) x$waic$waic),
    p_W = sapply(ml$cont_PROV$res,
                 function(x) x$waic$p.eff),
    marginal_loglik = sapply(ml$cont_PROV$res,
                      function(x) x$mlik[1, 1]),
    LPML = sapply(ml$cont_PROV$res,
                      function(x) sum(log(x$cpo$cpo), na.rm = TRUE))
  )

kable(
  contPROV_mcdat,
  digits = 0
)

ggplot(contPROV_mcdat[-1, ], aes(DIC, marginal_loglik)) +
  geom_point() + 
  geom_text(aes(label = rownames(contPROV_mcdat)[-1]), hjust = -.1)

PROVdat <- local({

  ## Predictions of the categorical model

  cat <- with(
    ml$cont_PROV$res$cat,
    rbind(summary.fixed[1:3, c('mean', '0.025quant', '0.975quant')],
          summary.fixed[4, 'mean'] +
            summary.fixed[1:3, c('mean', '0.025quant', '0.975quant')])
  )


  catxyear <- with(
    ml$cont_PROV$res$catxyear,
        summary.fixed[t(matrix(1:6, 2)),
                      c('mean', '0.025quant', '0.975quant')]
  )


  center.ranef <-
    sapply(
    1:2,
    function(x) 
      mean(ml$cont_PROV$spderes$cat$summary.values$mean[st.index$theta.group == x],
           na.rm = TRUE)) +
    mean(ml$cont_PROV$res$cat$summary.random$recode.cont$mean[gen_str$recoding],
           na.rm = TRUE)



  data.frame(
    Year = factor(rep(2012:2013, each = 3)),
    centering = rep(center.ranef, each = 3),
    PROV   = factor(rep(levels(ml$dat$PROV), 2)),
    cat = cat,
    cy = catxyear ) %>% 
    gather(model, y, cat.mean:cy.0.975quant) %>% 
    separate(model, c('model', 'q'), extra = 'merge') %>% 
    mutate(y = y + centering) %>% 
    spread(q, y)
})

ggplot(PROVdat, aes(PROV, mean, colour = model)) +
  geom_pointrange(aes(ymin = `0.025quant`, ymax = `0.975quant`),
                  position = position_dodge(width = 0.5)) +
  facet_wrap(~ Year)

## Back-transform the posterior means, into the prob. scale
## Note this is not the posterior mean in that scale, but approx.


tdat <- filter(mdat,
               variable %in% c('CD', 'CL')
               ,Year %in% c('2012', '2013')
) %>%
  spread(variable, value) %>% 
  transform(CLcont = factor(CL > 0, labels = c('0', '>0'))) %>%
  filter(!is.na(CL), !is.na(PROV)) %>% 
  droplevels() %>% 
  group_by(Year, PROV) %>% 
  summarise(y = mean(CL[CL > 0]),
            n = n())

ggplot(mutate(PROVdat, y = exp(mean)), aes(PROV, y)) +
  geom_bar(aes(fill = model), stat = 'identity', position = 'dodge') + 
  geom_point(dat = tdat, stat = 'identity') +
  facet_wrap(~ Year) +
  ylab("CL | CL > 0")

The provenance is clearly irrelevant.

Role of the Bud-flush

Here we test the relevance of the Bud-flush in explaining the Binary component.

## Run models for comparison
source('CL_cont_vs_BF_models.R')

ml$cont_BF$res <- c(lapply(ml$cont_BF$res, inla.qget))

## Summary of spatial results
ml$cont_BF$spderes <- lapply(ml$cont_BF$res[-1], inla.spde2.result, 'theta.cont', sp_str$spde)

# saveRDS(ml$cont_BF, file.path(path, 'cache', 'CL.ml.cont_BF.rds'))
ml$cont_BF <- RisingAshes:::BFmodels_CLcont

lapply(ml$cont_BF$res, summary)

BF_mcdat <- 
    data.frame(
    DIC = sapply(ml$cont_BF$res,
                 function(x) x$dic$dic),
    WAIC = sapply(ml$cont_BF$res,
                  function(x) x$waic$waic),
    marginal_loglik = sapply(ml$cont_BF$res,
                      function(x) x$mlik[1, 1]),
    LPML = sapply(ml$cont_BF$res,
                      function(x) sum(log(x$cpo$cpo), na.rm = TRUE))
  )

kable(
  BF_mcdat,
  digits = 0
)

ggplot(BF_mcdat[-1, ], aes(DIC, marginal_loglik)) +
  geom_point() + 
  geom_text(aes(label = rownames(BF_mcdat)[-1]), hjust = -.1)

BFdat <- local({

  nlevels <- 5
  ## Predictions of the categorical model

  cat <- with(
    ml$cont_BF$res$cat,
    rbind(summary.fixed[1:nlevels, c('mean', '0.025quant', '0.975quant')],
          summary.fixed[nlevels+1, 'mean'] +
            summary.fixed[1:nlevels, c('mean', '0.025quant', '0.975quant')])
  )


  catxyear <- with(
    ml$cont_BF$res$catxyear,
        summary.fixed[t(matrix(1:(2*nlevels), 2)),
                      c('mean', '0.025quant', '0.975quant')]
  )


  center.ranef <-
    sapply(
    1:2,
    function(x) 
      mean(ml$cont_BF$spderes$cat$summary.values$mean[st.index$theta.group == x],
           na.rm = TRUE)) +
    mean(ml$cont_BF$res$cat$summary.random$recode.cont$mean[gen_str$recoding],
           na.rm = TRUE)



  data.frame(
    Year = factor(rep(2012:2013, each = nlevels)),
    centering = rep(center.ranef, each = nlevels),
    BF   = factor(rep(levels(ml$dat$BF98), 2)),
    cat = cat,
    cy = catxyear ) %>% 
    gather(model, y, cat.mean:cy.0.975quant) %>% 
    separate(model, c('model', 'q'), extra = 'merge') %>% 
    mutate(y = y + centering) %>% 
    spread(q, y)
})

ggplot(BFdat, aes(BF, mean, colour = model)) +
  geom_pointrange(aes(ymin = `0.025quant`, ymax = `0.975quant`),
                  position = position_dodge(width = 0.5)) +
  facet_wrap(~ Year)

## Back-transform the posterior means, into the prob. scale
## Note this is not the posterior mean in that scale, but approx.

tdat <- filter(mdat,
               variable %in% c('CD', 'CL')
               ,Year %in% c('2012', '2013')
) %>%
  spread(variable, value) %>% 
  transform(CLcont = factor(CL > 0, labels = c('0', '>0'))) %>%
  filter(!is.na(CL), !is.na(BF98)) %>% 
  droplevels() %>% 
  group_by(Year, BF98) %>% 
  summarise(y = mean(CL[CL>0])) %>% 
  mutate(BF = BF98)

ggplot(mutate(BFdat,
              y = exp(mean)), aes(BF, y)) +
  geom_bar(aes(fill = model), stat = 'identity', position = 'dodge') + 
  geom_point(dat = tdat, stat = 'identity') +
  facet_wrap(~ Year) +
  ylab("CL | CL > 0")

Absolutely irrelevant.

Model selection continuous component

We have analysed the effects of the Basal Circumference (BC), the Provenance (PROV) and the Bud-flushing (BF) separately, on the continuous variable CL | CL > 0.

In all three cases, there is no evidence of a relationship between these variables and the response. We are not including them in the continous component of the model.

Mixture spatio-temporal model

For a measurement of CL $y_{ijk}$ taken in year $i$, at the location $j$ for the individual $k$, we assume that \begin{equation} \begin{aligned} \text{Pr}[y_{ijk} = 0]& = p_{ijk}, && \quad 0 < p_{ijk} < 1 \ \pi(y_{ijk} | y_{ijk} > 0)& = \text{Ga}(a_{ijk}, b_{ijk}), && \quad a_{ijk}, b_{ijk} > 0. \end{aligned} \end{equation} Or, equivalently, that the probability density of the response variable is \begin{equation} \pi(y | p, \alpha, \beta) = p \delta_0 + (1-p) \text{Ga}(y | a, b) I_{[y > 0]}, \end{equation} where $\delta_0$ is the Dirac delta function

This can be seen as a semi-continuous model for the observed data, built as a mixture of a Bernoulli distribution which gives the pattern of zeros and a Gamma distribution for the positive values. We define a hierarchical model for the parameters $p$, $a$ and $b$ using appropriate link functions of the expected values of the respective distributions. Specifically, calling $\mu = E(y | y > 0) = \frac{a}{b}$, we define two linear predictors \begin{equation} \begin{split} \text{logit}(p_{ijk}) & = \text{X}{ik}^{(1)} \beta^{(1)} + s{ij}^{(1)} + a_k^{(1)}\ \text{log}(\mu_{ijk}) & = \text{X}{ik}^{(2)} \beta^{(2)}+ s{ij}^{(2)} + a_k^{(2)}. \end{split} \end{equation} where $\text{X}{ij}^{(\cdot)}$ are matrices with explanatory covariates specific for each linear predictor potentially depending on the year and the individual; $s{ij}$ is a structured spatio-temporal random effect and $a_k$ is an structured additive-genetic random effect at individual level.

The spatio-temporal structure is built as the Kronecker product of separate temporal and spatial zero-mean Gaussian processes. The temporal process is simply defined by a correlation parameter $\rho_t$ between consecutive observations, while the spatial process is defined by a Matérn covariance function with parameters of shape (or smoothness) $\nu$, spatial scale $\kappa$ and precision (i.e. inverse variance) $\tau$. The smoothness parameter is fixed to $\nu = 1$ for convenience. The spatial scale parameter is associated with the effective range of the spatial process, so that the correlation between locations at a distance $\rho_s = \sqrt{8}/\kappa$ is approximately $0.13$. Finally, the marginal variance of the spatial process is given by $\sigma_s^2 = (4 \pi\kappa^2\tau^2)^{-1}$.

Finally, the structured additive-genetic effect follows a zero-mean multivariate Normal distribution with a known covariance structure given by the family kinship. Specifically, the covariance matrix is $G = \sigma_a^2 A$, where $\sigma_a^2$ is the additive genetic variance in the base population and the additive-genetic structure matrix $A$ has elements $A_{kk'} = 2 \Theta_{kk'}$; i.e., twice the coefficient of coancestry between the individuals $k$ and $k'$ (see for example, Lynch and Walsh 1998, p. 756).

The results from this mixture model are equivalent to those obtained by running two separate models, one for each component. The difference is in the evaluation of the likelihood and other measures like DIC and WAIC. Running both models jointly allows for model comparison with other settings.

## Dataset
ml$mst$dat <- 
  ml$dat[, c('Seq', 'FAM', 'PROV', 'BF98', 'Year', 'BC')] %>% 
  transform(
    Year2013 = as.numeric(Year == 2013),
    BC = cut(BC, c(0, 30, 45, 60, 75, 90, 150), ordered = TRUE),
    BF_2 = as.numeric(BF98 == 2),
    BF_3 = as.numeric(BF98 == 3),
    BF_4 = as.numeric(BF98 == 4),
    BF_5 = as.numeric(BF98 == 5),
    recode = gen_str$recoding,
    FAM = as.numeric(FAM),
    CL0 = as.numeric(ml$dat$CL == 0))

ml$mst$dat$CL_gt0           <- ml$dat$CL
ml$mst$dat$CL_gt0[ml$dat$CL == 0] <- NA

ml$mst$dat$yt <- ml$mst$dat$CL_gt0
ml$mst$dat$yt[ml$dat$CL == 1] <- .98   # forbiden value for a beta


## Binomial stack
stk.bin <-
  inla.stack(
    data    = list(y = cbind(as.numeric(ml$dat$CL == 0), NA)),
    A       = list(st.Amat, 1, st.Amat, 1),
    effects = list(suffix(st.index, 'bin'),
                   suffix(ml$mst$dat, 'bin'),
                   suffix(st.index, 'cont', makeNA = TRUE),
                   suffix(ml$mst$dat, 'cont', makeNA = TRUE)),
    tag     = 'bin')


## continuous stack

stk.cont <-
  inla.stack(data    = list(y = cbind(NA, ml$mst$dat$yt)),
             A       = list(st.Amat, 1, st.Amat, 1),
             effects = list(suffix(st.index, 'bin', makeNA = TRUE),
                            suffix(ml$mst$dat, 'bin', makeNA = TRUE),
                            suffix(st.index, 'cont'),
                            suffix(ml$mst$dat, 'cont')),
             tag     = 'cont')

ml$mst$stk <- inla.stack(stk.bin, stk.cont)


## Prior on the precision of the overdispersion term
## We use a Gamma with shape 0.5 and rate 0.2, which gives a
## prior expectation (in the log scale) of 2.5 and a reasonably ample
## variance of 12.5.
## Furthermore, we checked that sensitivity to this prior specification is very
## limited.
prec.od <- list(fixed = FALSE,
                initial = 1,
                param = c(.5, .2))


## Formula
ml$mst$fml <- y ~ 0 +
  Year.bin +
  # Year.bin:BC.bin +BF_2.bin + BF_3.bin + BF_4.bin + BF_5.bin +
  f(theta.bin, model=sp_str$spde,
    group = theta.group.bin,  # variable in the mesh space
    control.group = list(model = 'exchangeable')) +
  f(recode.bin, model = "generic0", hyper = list(theta = prec.A),
    constr = TRUE, Cmatrix = Cmat) +
  # f(Seq.bin, model = "iid", hyper = list(theta = prec.od)) + # only in the continous part!!
  Year.cont +
  f(theta.cont, model=sp_str$spde,
    group = theta.group.cont,  # variable in the mesh space
    control.group = list(model = 'exchangeable')) +
  f(recode.cont, model = "generic0", hyper = list(theta = prec.A),
    constr = TRUE, Cmatrix = Cmat) 
# + f(Seq.cont, model = "iid", hyper = list(theta = prec.A))  # Causes overfitting


# Use the right link function for the prediction when the response is NA
link.NA <- rep(NA, nrow(ml$dat)*2L)
nas.idx.bin <- which(is.na(ml$dat$CL))
nas.idx.cont <- which(is.na(ml$mst$dat$yt))   # there were additional NAs
link.NA[nas.idx.bin] <- 1
link.NA[nrow(ml$dat) + nas.idx.cont] <- 2

## Introduce the "residual" variation into a latent 
## overdispersion term and fix the likelihood precision
## to a high value
ml$mst$res <- inla(
  formula = ml$mst$fml,
  data = inla.stack.data(ml$mst$stk),
  family = list("binomial", "gamma"),
  Ntrials = 1,
  control.family = list(
    list(),
    list(
      hyper = list(theta = prec.od)
      # hyper = list(theta = list(fixed = TRUE, initial = 10))
      # causes overfitting
    )
  ),
  ## Correction of the Laplace Approximation
  ## for Binomial data with many zeros
  ## http://arxiv.org/abs/1503.07307
  control.inla = list(
    correct = TRUE,
    correct.factor = 10),
  control.predictor = list(
    A = inla.stack.A(ml$mst$stk),
    compute = TRUE, link = link.NA),
  control.fixed = list(
    expand.factor.strategy = 'inla'), # account for NAs in factor
  control.compute = list(
    dic = TRUE,
    waic = TRUE,
    cpo = TRUE,
    po = TRUE,
    config = TRUE)
  #, control.family = list(hyper = list(prec = prec.A))
  # , inla.call = 'submit'
)

ml$mst$res <- inla.qget(ml$mst$res)

ml$mst$res <- RisingAshes:::mst_CL

## Summary of spatial results
ml$mst$spde.bin <- 
  inla.spde2.result(ml$mst$res, 'theta.bin', sp_str$spde)
ml$mst$spde.cont <-
  inla.spde2.result(ml$mst$res, 'theta.cont', sp_str$spde)

Model summary

summary(ml$mst$res)

R-INLA package version r ml$mst$res$version$R.INLA

Fitted vs. Observations

Note that the fact that the prediction is not centered around the observed values is a consequence of the discrete nature of the response variable. It would be unlikely to predict a value under 1 or above 3.16, when these are the extreme observed values.

bin.idx <- inla.stack.index(ml$mst$stk, tag = 'bin')$data
cont.idx <- inla.stack.index(ml$mst$stk, tag = 'cont')$data

## If I didn't use link = 2 in the INLA call:
fitted.bin <- ml$mst$res$summary.fitted.values$mean[bin.idx]
for( i in nas.idx.bin ) {
  fitted.bin[i] <- inla.emarginal(inv.logit, ml$mst$res$marginals.fitted.values[[bin.idx[i]]])
}
fitted.cont <- ml$mst$res$summary.fitted.values$mean[cont.idx]
for( i in nas.idx.cont ) {
  fitted.cont[i] <- inla.emarginal(inv.logit, ml$mst$res$marginals.fitted.values[[cont.idx[i]]])
}

fvso.dat <- rbind(
  data.frame(
    component = 'bin',
    Observed = as.numeric(ml$dat$CL == 0),
    Predicted = fitted.bin,
    Residuals = as.numeric(ml$dat$CL == 0) - fitted.bin,
    Year = ml$dat$Year
  ),
  data.frame(
    component = 'cont',
    Observed = ml$mst$dat$yt,
    Predicted = fitted.cont,
    Residuals = ml$mst$dat$yt - fitted.cont,
    Year = ml$dat$Year
  )
)

g <- ggplot(fvso.dat, aes(Observed, Predicted)) +
  geom_violin(aes(group = Observed),
              data = subset(fvso.dat, component == 'bin'),
              fill = 'grey80') + 
  geom_point(aes(x=Observed),
             data = subset(fvso.dat, component == 'cont')) + 
  geom_abline(intercept = 0, slope = 1, col = 'darkgray',
              data = subset(fvso.dat, component == 'cont')) +
  facet_wrap(~ component, scales = 'free_x')

## Tweak the x-axis for the left facet
## ggplot under the hood:
## http://zevross.com/blog/2014/11/20/under-the-hood-of-ggplot2-graphics-in-r/
gTable <- ggplot_gtable(ggplot_build(g))
## Determine the right grob
# pdf("grobs.pdf")
# for(i in 1:length(gTable$grobs)){
#   grid.draw(gTable$grobs[[i]])
#   grid.text(i, x = unit(0.1, "npc"), y = unit(0.1, "npc"))
#   grid.newpage()
# }
# dev.off()
gTable$grobs[[8]]$children[2]$axis$grobs[[2]]$children[[1]]$label <- 
  c('', '> 0', '', 'Zero', '')
gTable$grobs[[8]]$children[2]$axis$grobs[[1]]$y <- 
  unit(0, units = 'cm')
grid.draw(gTable)

ggplot(cbind(fvso.dat, BC = ml$dat$BC), aes(BC, Residuals, colour = Year)) +
  geom_point() +
  facet_wrap(~ component)

BC = cut(ml$dat$BC, c(0, 30, 45, 60, 75, 90, 150), ordered = TRUE)

cbind(fvso.dat, BC) %>% 
  group_by(BC, Year, component) %>% 
  summarise(mresid = mean(Residuals, na.rm = TRUE),
            N = n()) %>% 
  ggplot(aes(BC, mresid, colour = Year)) +
  geom_point(aes(size = N)) +
  facet_wrap(~ component) +
  ylab("Mean residuals")

It is certainly not evident at all from this diagnostics plot that there is a missing relationship between CL and BC. The strong differences for circumferences under 30 cm are based on very few points, and might well have arosen by chance.

ggplot(cbind(fvso.dat, BC = ml$dat$PROV),
       aes(BC, Residuals, colour = Year)) +
  geom_violin() +
  facet_wrap(~ component)

Nor for the case of the provenances.

ggplot(cbind(fvso.dat, BC = ml$dat$BF98), aes(BC, Residuals, colour = Year)) +
  geom_violin() +
  facet_wrap(~ component)

Nor for the flushing.

Genetic effects

Here we plot the posterior distribution for the individual breeding values. Family effect included implicitly.

Not very useful, isn't it?

hpds.bin <- sapply(ml$mst$res$marginals.random$recode.bin,
                   function(x) inla.hpdmarginal(.95, x))

ml$mst$blups.bin <- data.frame(
  Seq = dat$Seq,
  Family = dat$FAM,
  mean = ml$mst$res$summary.random$recode.bin$mean[gen_str$recoding],
  ymin = hpds.bin[1, gen_str$recoding],
  ymax = hpds.bin[2, gen_str$recoding])

ggplot(ml$mst$blups, aes(x = mean, y = mean, ymin = ymin, ymax = ymax)) + 
  geom_pointrange() + 
  facet_wrap(~ Family) +
  coord_flip() +
  labs(x=NULL, y=NULL)

Posterior distributions of spatial field's characteristics

The spatial autocorrelation extends over a range of about 20 rows/columns, with a significant variability.

stopifnot(exists('spde.bin', ml$mst))
res.spde.bin.marginals <- 
  ldply(tail(ml$mst$spde.bin, 2),
        function(x) data.frame(x[[1]]))
res.spde.bin.marginals <- 
  transform(res.spde.bin.marginals,
            .id=factor(res.spde.bin.marginals$.id,
                       labels = c('Range', 'Variance')))
res.spde.bin.marginals <- 
  rbind(
    cbind(type = 'Posterior',
          res.spde.bin.marginals),
    cbind(type = 'Prior', .id = 'Range',
          as.data.frame(sp_str$priors$rho)),
    cbind(type = 'Prior', .id = 'Variance',
          as.data.frame(sp_str$priors$sigma))
  )

# qplot(x, y, data=res.spde.bin.marginals, geom='line', col = type) +
#   facet_wrap(~.id, scales='free')

stopifnot(exists('spde.cont', ml$mst))
res.spde.cont.marginals <- 
  ldply(tail(ml$mst$spde.cont, 2),
        function(x) data.frame(x[[1]]))
res.spde.cont.marginals <- 
  transform(res.spde.cont.marginals,
            .id=factor(res.spde.cont.marginals$.id,
                       labels = c('Range', 'Variance')))
res.spde.cont.marginals <- 
  rbind(
    cbind(type = 'Posterior',
          res.spde.cont.marginals),
    cbind(type = 'Prior', .id = 'Range',
          as.data.frame(sp_str$priors$rho)),
    cbind(type = 'Prior', .id = 'Variance',
          as.data.frame(sp_str$priors$sigma)[sp_str$priors$sigma$x < 1, ]))
# qplot(x, y, data=res.spde.cont.marginals, geom='line', col = type) +
#   facet_wrap(~.id, scales='free')

hyper.dat <- 
  rbind(
    cbind(
      component = 'bin',
      res.spde.bin.marginals),
    cbind(component = 'cont',
          res.spde.cont.marginals)
  )

grid.arrange(
  ggplot(filter(hyper.dat, .id == 'Range'), aes(x, y, lty = type)) + 
    geom_line(show.legend = FALSE) + 
    facet_wrap(component ~ .id, scales = 'free_x', ncol = 1) +
    labs(x=NULL, y=NULL) +
    xlim(0, 150) +
    theme(legend.title = element_blank()),
  ggplot(filter(hyper.dat, .id == 'Variance'), aes(x, y, lty = type)) + 
    geom_line() + 
    facet_wrap(component ~ .id, scales = 'free_x', ncol = 1) +
    labs(x=NULL, y="") +
    # xlim(0, 150) +
    theme(legend.title = element_blank()),
  nrow = 1,
  widths = 3:4
)

Posterior spatiotemporal field

The surfaces for the years 2012 and 2013 are remarkably similar in both components. We should try a model with one single spatial effect for both years.

This indicates, unlike the Crown Dieback variable, that there is no spatial spread of this symptom.

proj.vec = inla.mesh.projector(sp_str$mesh, loc = sp_str$loc)
sp.f <- function(x) inla.mesh.project(proj.vec, field = x)

ml$mst$post_spat.bin <- ldply(1:2, function(x) 
  data.frame(X = dat$X,
             Y = dat$Y,
             Year = levels(ml$dat$Year)[x],
             PostMean= sp.f(ml$mst$spde.bin$summary.values$mean[st.index$theta.group == x])))


sp.bin <- ggplot(ml$mst$post_spat.bin, aes(X, Y)) +
  geom_tile(aes(fill = PostMean, color = PostMean)) +
  coord_fixed() +
  labs(x='', y='Field row') +
  facet_wrap(~ Year) +
  scale_fill_gradient2(low='#832424FF', high='#3A3A98FF', name = 'Binary\npost.\nmean', space = "Lab") +
  scale_color_gradient2(low='#832424FF', high='#3A3A98FF', name = 'Binary\npost.\nmean', space = "Lab")

proj.vec = inla.mesh.projector(sp_str$mesh, loc = sp_str$loc)
sp.f <- function(x) inla.mesh.project(proj.vec, field = x)

ml$mst$post_spat.cont <- ldply(1:2, function(x) 
  data.frame(X = dat$X,
             Y = dat$Y,
             Year = levels(ml$dat$Year)[x],
             PostMean= sp.f(ml$mst$spde.cont$summary.values$mean[st.index$theta.group == x])))


sp.cont <- ggplot(ml$mst$post_spat.cont, aes(X, Y)) +
  geom_tile(aes(fill = PostMean, color = PostMean)) +
  coord_fixed() +
  labs(x='Field column', y='Field row') +
  facet_wrap(~ Year) +
  scale_fill_gradient2(low='#3A3A98FF', high='#832424FF', name = 'Gamma\npost.\nmean', space = "Lab") +
  scale_color_gradient2(low='#3A3A98FF', high='#832424FF', name = 'Gamma\npost.\nmean', space = "Lab")

grid.arrange(sp.bin, sp.cont)

## shades of gray version for the journal
comp_plot <- 
  grid.arrange(
    ggplot(ml$mst$post_spat.bin, aes(X, Y)) +
      geom_tile(aes(fill = PostMean, color = PostMean)) +
      coord_fixed() +
      labs(x='', y='Field row') +
      facet_wrap(~ Year) +
      scale_fill_gradient(
        low = 'black', high = 'white',
        name = 'Binary\npost.\nmean', space = "Lab") +
      scale_color_gradient(
        low = 'black', high = 'white',
        name = 'Binary\npost.\nmean', space = "Lab")
    ,
    ggplot(ml$mst$post_spat.cont, aes(X, Y)) +
      geom_tile(aes(fill = PostMean, color = PostMean)) +
      coord_fixed() +
      labs(x='Field column', y='Field row') +
      facet_wrap(~ Year) +
      scale_fill_gradient(
        low = 'white', high = 'black',
        name = 'Gamma\npost.\nmean', space = "Lab") +
      scale_color_gradient(
        low = 'white', high = 'black',
        name = 'Gamma\npost.\nmean', space = "Lab")
  )

## A list of 942 * 2 * 2 = 3768 marginals
## corresponding for the SPDE vertices of
## bin-2012; bin-2013; cont-2012; cont-2013
## in that specific order
spmar <- c(ml$mst$spde.bin$marginals.values,
           ml$mst$spde.cont$marginals.values)

## mean values of linear predictor for each Year
latent_means <- ml$mst$res$summary.fixed[, 'mean']

## ST marginals centered on the Year's mean
spmar.shifted <- mapply(
  shift_marginal, 
  spmar, 
  latent_means[c(st.index$theta.group,
                 2 + st.index$theta.group)],
  SIMPLIFY = FALSE)

## ST means on the original scale for each Year and component
spmeans.oscale <- c(
  sapply(
    spmar.shifted[seq_len(2*sp_str$mesh$n)],             # bin marginals
    function(x) inla.emarginal(inv.logit, x)
  ),
  sapply(
    spmar.shifted[2*sp_str$mesh$n + seq_len(2*sp_str$mesh$n)],  # cont marginals
    function(x) inla.emarginal(exp, x)
  )
)

ml$mst$post_spat_oscale <- ldply(
  1:2, function(x) 
    data.frame(
      X = dat$X,
      Y = dat$Y,
      Year = (2012:2013)[x],
      PostMean_bin= sp.f(spmeans.oscale[seq_len(2*sp_str$mesh$n)][st.index$theta.group == x]),
      PostMean_cont= sp.f(spmeans.oscale[2*sp_str$mesh$n+seq_len(2*sp_str$mesh$n)][st.index$theta.group == x])
    )
)


sp.bin <- ggplot(ml$mst$post_spat_oscale, aes(X, Y)) +
  geom_tile(aes(fill = PostMean_bin, color = PostMean_bin)) +
  coord_fixed() +
  labs(x='', y='Field row') +
  facet_wrap(~ Year) +
  scale_fill_viridis(name = 'Binary\npost.\nmean') +
  scale_color_viridis(name = 'Binary\npost.\nmean') 


sp.cont <- ggplot(ml$mst$post_spat_oscale, aes(X, Y)) +
  geom_tile(aes(fill = PostMean_cont, color = PostMean_cont)) +
  coord_fixed() +
  labs(x='Field column', y='Field row') +
  facet_wrap(~ Year) +
  scale_fill_viridis(name = 'Gamma\npost.\nmean') +
  scale_color_viridis(name = 'Gamma\npost.\nmean')

In the following figure we can see the mean posterior ST effect of both model components, by Year, in the original scale.

grid.arrange(sp.bin, sp.cont)

Phenotype vs. Predicted Breeding Values (PBV)

plotdat <- ml$dat %>%
  dplyr::select(Seq, Year, CL) %>%
  spread(Year, CL) %>%
  inner_join(ml$mst$blups.bin, by = 'Seq') %>%
  mutate(Infection = factor((`2012` > 0) + 2*(`2013` > 0),
                            labels = c('Not infec.',
                                       'Only 2012',
                                       'Only 2013',
                                       'Both')))

ggplot(plotdat, aes(x = mean, fill = Infection)) +
  geom_histogram(binwidth = diff(range(plotdat$mean))/30) +
  scale_fill_grey(na.value = 'white') +
  xlab('Predicted Breeding Value')

## Simplified version for journal
ggplot(plotdat[!is.na(plotdat$Infection) & plotdat$Infection != 'Only 2012', ], aes(x = mean)) +
  geom_histogram(binwidth = diff(range(plotdat$mean))/30) +
  facet_grid(Infection~.) +
  xlab('Predicted Breeding Value')

Heritability

In principle, we have two heritabilities here, corresponding to two different traits: the first explains the pattern of the zeros, and the second explains the degree of infection, given that there is infection. Can be thought of like resistance and sensitivity.

Both with a Binomial and Gamma likelihoods, the phenotypic variance is a function of the mean, which means that it varies from observation to observation. Although the additive-genetic variance(s) is common to the population, the phenotypic variance can not be decomposed additively into genetic and residual components. In this context some funny things can take place. For example, a decrease in the phenotypic variance as a consequence of an increase in the genetic variance.

Heritability in the latent/observed scale

In order to compute the posterior density of the heritability, we sample from the posterior distributions of the hyperparameters of the model, and build a posterior sample of heritabilities with the formula [@Nakagawa10, @Vazquez09] $$ h^2 = \frac{\sigma^2_a}{\sigma^2_a + \omega \pi^2/3}, $$ where $\omega$ equals $1$ in the case of the Binary trait.

For the Gamma component, we followed the general approach described in [@Villemereuil15].

N.samples <- 5e3

hyper.dat <- rbind(
  cbind(component = 'bin',
        res.spde.bin.marginals),
  cbind(component = 'cont',
        res.spde.cont.marginals)
)


sbin.mar <- filter(hyper.dat,
                   component == 'bin',
                   type == 'Posterior',
                   .id == 'Variance')
sbet.mar <- filter(hyper.dat,
                   component == 'cont',
                   type == 'Posterior',
                   .id == 'Variance')
hyperpar.sample <- inla.hyperpar.sample(N.samples, ml$mst$res)

var.sample <- data.frame(
  sbin = inla.rmarginal(N.samples, sbin.mar[, c('x', 'y')]),
  sbet = inla.rmarginal(N.samples, sbet.mar[, c('x', 'y')]),
  gbin  = 1/hyperpar.sample[, grep('recode.bin',
                                   colnames(hyperpar.sample))],
  gbet  = 1/hyperpar.sample[, grep('recode.cont',
                                   colnames(hyperpar.sample))],
  omega = 1/(1 + hyperpar.sample[, grep('Gamma observations',
                                        colnames(hyperpar.sample))])
)


## Heritability in the observed scale
## Approach following Villemereuil et al. 2015

# Samples from fitted values (i.e. exp(latent))

## Numerical problems with some very spliked marginals
## for the binary component. E.g.:
# tm <- ml$mst$res$marginals.fitted.values[inla.stack.index(ml$mst$stk,
#                                                       'bin')$data][[31]]
# plot(tm, type = 'l')
# inla.qmarginal(.1, tm)
# str(inla.smarginal(tm))
# plot(inla.smarginal(tm), type = 'l')

# The following discussion can be useful to overcome these numerical problems:
# https://groups.google.com/forum/embed/?parenturl=http%3A%2F%2Fwww.r-inla.org%2Fcomments-1&service=jotspot&ul=1&theme=default&place=forum%2Fr-inla-discussion-group&showsearch=true#!topic/r-inla-discussion-group/9B8CazRaik4)

## Other marginals give samples out of scale. E.g.: #16
# tm <- ml$mst$res$marginals.fitted.values[inla.stack.index(ml$mst$stk,
#                                                       'bin')$data][[31]]
# stopifnot(inla.rmarginal(N.samples, tm) > 0)

## Return a sample of N values from the posterior
## marginals of the fitted.values of the given component
sample_fit <- function(comp, mar, stk, N) {
  sapply(mar[inla.stack.index(stk, comp)$data],
         function(x) {
           tryCatch(inla.rmarginal(N, x),
                    error = function(e) 
                      rep(NA, N))
           }
  )
}

## Sample the fitted values for each component
fit.sample <- lapply(
  c(bin = 'bin', cont = 'cont'),
  sample_fit,
  mar = ml$mst$res$marginals.fitted.values,
  stk = ml$mst$stk,
  N   = N.samples
)

# Due to numerical problems in some marginals, I replace the 
# random sampling by the mode of the posterior
idx <- lapply(fit.sample, apply, 2, function(x) all(is.na(x)))
stopifnot(sapply(idx, sum) < sapply(idx, length) / .05) # < 5 %
fit.sample$bin[, idx$bin] <- 
  rep(
    ml$mst$res$summary.fitted.values$mean[
      inla.stack.index(ml$mst$stk, 'bin')$data][
        idx$bin],
    each = N.samples
  )

# Auxiliar diagnostic function
plot_sample <- function(x) {
  x %>% 
    as.data.frame() %>% 
    gather(component) %>% 
    ggplot(aes(value)) +
    geom_histogram() +
    facet_wrap(~component, scales = 'free')
}

# Average derivative of expected value wrt latent value (Eqs. 16-17)
Psi.sample <- lapply(fit.sample, apply, 1, mean)
# plot_sample(Psi.sample)

# Expected-scale phenotypic variance (Eq. 7)
V_P_exp.sample <- lapply(fit.sample, apply, 1, var)
# plot_sample(V_P_exp.sample)

# Distribution-specific variance:
# (MC integrated wrt latent)
# Gamma distribution: mu^2/phi
# Bernoulli distribution: p(1-p)
phi.sample <- hyperpar.sample[, grep('Gamma observations',
                                     colnames(hyperpar.sample))]
V_P_spec.sample <- list(
  bin = apply(fit.sample$bin*(1-fit.sample$bin), 1, mean),
  cont = apply(fit.sample$cont^2/phi.sample, 1, mean)
)
# plot_sample(V_P_spec.sample)

# Observed-scale phenotypic variance
V_P_obs <- mapply(`+`, V_P_exp.sample, V_P_spec.sample,
                  SIMPLIFY = FALSE)
# plot_sample(V_P_obs)

# Additive-genetic variance (latent scale)
V_A <- list(
  bin = 1/hyperpar.sample[, grep('recode.bin', colnames(hyperpar.sample))],
  cont = 1/hyperpar.sample[, grep('recode.cont', colnames(hyperpar.sample))]
)
# plot_sample(V_A)

# Additive-genetic variance (observed/expected phenotypic scales)
V_A_obs <- mapply(function(x, y) x^2*y,
                  Psi.sample,V_A,
                  SIMPLIFY = FALSE)
# plot_sample(V_A_obs)

## Method a: dividing by ST variance
## Method b: excluding ST variance
ml$mst <- c(ml$mst[-grep('h2', names(ml$mst))],
            transmute(var.sample,
                      h2.sample.bin.b = gbin/(gbin + pi^2/3),
                      # h2.sample.cont.b = gbet/(gbet + omega*pi^2/3),
                      # This formula applied for beta, but not for gamma
                      h2.sample.bin.a = gbin/(gbin + sbin + pi^2/3),
                      # h2.sample.bin.a = V_A_obs$bin/V_P_obs$bin,
                      # This give very high values. Numerical issues?
                      h2.sample.cont.a = V_A_obs$cont/V_P_obs$cont))


plotdat <- as.data.frame(ml$mst[grep('h2', names(ml$mst))]) %>% 
  # dplyr::rename(bin = h2.sample.bin, cont = h2.sample.cont) %>% 
  gather(component, x) %>% 
  mutate(component = factor(gsub('h2.sample.', '', component))) %>% 
  separate(component, c('component', 'spatial_variance'), convert = TRUE)


ggplot(plotdat, aes(x, colour = component)) +
  geom_density() +
  facet_grid(`spatial_variance`~.)

I don't have a good way of computing the heritability for the Gamma component excluding the ST variance. On the other hand, the approach by [@Villemereuil15] gave numerical issues for the binary compoent.

(plotdat %>% 
  group_by(component, spatial_variance) %>% 
  dplyr::summarise(mean = mean(x),
                   ymin = quantile(x, probs = c(.025)),
                   ymax = quantile(x, probs = c(.975))) -> plotdatsum)

ggplot(plotdatsum, aes(component, mean, ymin = ymin, ymax = ymax)) +
  geom_pointrange() +
  facet_grid(`spatial_variance`~.) +
  coord_flip() +
  ylim(c(0,1)) +
  ylab('Posterior heritability')

Conclusions

Here is a summary of the model comparison conducted separately for the binary and continuous components.

## Gathering of previous results
bin_mcdat <- rbind(
  ## Null model and 5 models for BC
  binBC_mcdat[-(c(1,4)), -6],
  ## 2 models for PROV
  binPROV_mcdat[3:4, -6]
)

## Model numbering
rownames(bin_mcdat) <- NULL
bin_mcdat$Model <- paste0('M', 
                          c(0, paste(1, 1:5, sep = '.'),
                            paste(2, 1:2, sep = '.')))

## Model descriptions
bin_mcdat$Specificity <-   
  c('Reference',
    'BC categories',
    'BC cat. x Year',
    'BC linear',
    'BC quadratic',
    'BC non-parametric',
    'PROV',
    'PROV x Year')

## Model components
bin_mcdat$`Linear predictor` <- 
  c("Year + _IBV_ + _ST_",
    "Year + BC + _IBV_ + _ST_",
    "Year * BC + _IBV_ + _ST_",
    "Year + $\\beta$ BC + _IBV_ + _ST_",
    "Year + $\\beta$ BC + $\\beta$ BC^2 + _IBV_ + _ST_",
    "Year + f(BC) + _IBV_ + _ST_",
    "Year + PROV + _IBV_ + _ST_",
    "Year * PROV + _IBV_ + _ST_"
  )

## -2 log-likelihood
bin_mcdat <- transform(bin_mcdat,
                       'Dev.' = -2*marginal_loglik)


## Reorder and select columns
bin_mcdat <- bin_mcdat[, c(6:8, 1:4, 9)]

kable(bin_mcdat, digits = 0,
      caption="Model selection criteria for the binary component")

## Gathering of previous results
cont_mcdat <- rbind(
  cont_comp[, -6],
  ## Null model and 5 models for BC
  contBC_mcdat[-(c(1,4)), -6],
  ## 2 models for PROV
  contPROV_mcdat[3:4, -6]
)

## Model numbering
rownames(cont_mcdat) <- NULL
cont_mcdat$Model <- paste0('M', 
                           c(paste(0, 1:3, sep = '.'), 
                             paste(1, 1:5, sep = '.'),
                             paste(2, 1:2, sep = '.')))

## Model descriptions
cont_mcdat$Specificity <-   
  c('Beta likelihood',
    'Gamma likelihood',
    'Reference',
    'BC categories',
    'BC cat. x Year',
    'BC linear',
    'BC quadratic',
    'BC non-parametric',
    'PROV',
    'PROV x Year')

## Model components
cont_mcdat$`Linear predictor` <- 
  c("BF98 + PROV + Year*BC",
    "BF98 + PROV + Year*BC",
    "Year + _IBV_ + _ST_",
    "Year + BC + _IBV_ + _ST_",
    "Year * BC + _IBV_ + _ST_",
    "Year + $\\beta$ BC + _IBV_ + _ST_",
    "Year + $\\beta$ BC + $\\beta$ BC^2 + _IBV_ + _ST_",
    "Year + f(BC) + _IBV_ + _ST_",
    "Year + PROV + _IBV_ + _ST_",
    "Year * PROV + _IBV_ + _ST_"
  )

## -2 log-likelihood
cont_mcdat <- transform(cont_mcdat,
                       'Dev.' = -2*marginal_loglik)


## Reorder and select columns
cont_mcdat <- cont_mcdat[, c(6:8, 1:4, 9)]

kable(cont_mcdat, digits = 0,
      caption="Model selection criteria for the continuous component")

References

famuvie/2016_RisingAshes documentation built on May 16, 2019, 10:02 a.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

famuvie/2016_RisingAshes
Companion to the Paper "Rising Out of the Ashes"

In famuvie/2016_RisingAshes: Companion to the Paper "Rising Out of the Ashes"

Models for the binary component

Role of the Basal Circumference

Role of the Provenance

Role of the Bud-flush

Model selection binary component

Models for the continuous component

Role of the Basal Circumference

Role of the Provenance

Role of the Bud-flush

Model selection continuous component

Mixture spatio-temporal model

Model summary

Fitted vs. Observations

Genetic effects

Posterior distributions of spatial field's characteristics

Posterior spatiotemporal field

Phenotype vs. Predicted Breeding Values (PBV)

Heritability

Heritability in the latent/observed scale

Conclusions

References

R Package Documentation

Browse R Packages

We want your feedback!

famuvie/2016_RisingAshes Companion to the Paper "Rising Out of the Ashes"

In famuvie/2016_RisingAshes: Companion to the Paper "Rising Out of the Ashes"

Models for the binary component

Role of the Basal Circumference

Role of the Provenance

Role of the Bud-flush

Model selection binary component

Models for the continuous component

Role of the Basal Circumference

Role of the Provenance

Role of the Bud-flush

Model selection continuous component

Mixture spatio-temporal model

Model summary

Fitted vs. Observations

Genetic effects

Posterior distributions of spatial field's characteristics

Posterior spatiotemporal field

Phenotype vs. Predicted Breeding Values (PBV)

Heritability

Heritability in the latent/observed scale

Conclusions

References

R Package Documentation

Browse R Packages

We want your feedback!

famuvie/2016_RisingAshes
Companion to the Paper "Rising Out of the Ashes"