• R-package implementing statistical models specifically suited for forest genetic resources analysts.
• Ultimately Mixed Models, but not necessarily easy to implement and use

• breedR acts as an interface which provides the means to:

  1. Combine any number of these models as components of a larger model
  2. Compute automatically incidence and covariance matrices from a few input parameters
  3. Fit the model
  4. Plot data and results, and perform model diagnostics

• Project web page http://famuvie.github.io/breedR/
  • download .zip or .tar.gz source files
  • install.packages('famuvie-breedR-*', type = 'source', repos = NULL)
• GitHub dev-site https://github.com/famuvie/breedR
  • if( !require(devtools) ) install.packages('devtools')
  • devtools::install_github('famuvie/breedR')
• CRAN not possible… ask I. Misztal
  • install.packages('breedR')

• Package's help: help(package = breedR)
  • Help pages ?remlf90
  • Code demos demo(topic, package = 'breedR') (omit topic for a list)
  • Vignettes vignette(package = 'breedR') (pkg and wiki)
• Wiki pages
  • Guides, tutorials, FAQ
• Mailing list http://groups.google.com/group/breedr
  • Questions and debates about usage and interface
• Issues page
  • Bug reports
  • Feature requests

• breedR is FOSS. Licensed GPL-3
  • RShowDoc('LICENSE', package = 'breedR')
• You can use and distribute breedR for any purpose
• You can modify it to suit your needs
  • we encourage to!
  • please consider contributing your improvements
  • you can distribute your modified version under the GPL
• However, breedR makes (intensive) use of the BLUPF90 suite of Fortran programs, which are for free but not free (remember CRAN?)

• Bayesian inference

• Multi-trait support

• Genotype\(\times\)Environment interaction *

• Support for longitudinal data

• Currently, only frequentist inference is supported via REML estimation of variance components.

• The function remlf90(), provides an interface to both REMLF90 and AIREMLF90 functions in the BLUPF90 suite of Fortran programs.

• Type ?remlf90 for details on the syntax

• It's on the roadmap by the end of this year

• Will use a gibbs sampler from BLUPF90, and possibly also INLA

• The interface will change a bit, separating the model specification from the fit

Specify the genetic group gg as an unstructured random effect using the standard formulas in R

\[ \begin{aligned} \mathrm{phe}_X = & \mu + Z \mathrm{gg} + \varepsilon \\ \mathrm{gg} \sim & N(0, \sigma_{\mathrm{gg}}^2) \\ \varepsilon \sim & N(0, \sigma_\varepsilon^2) \end{aligned} \]

res <- remlf90(fixed  = phe_X ~ 1,
               random = ~ gg,
               data   = globulus)

## No specification of initial variances.
##      Using default value of 1 for all variance components.
##      See ?breedR.getOption.

To avoid the notification, initial values for all the variance components must be made explicit using the argument var.ini:

res <- remlf90(fixed = phe_X ~ 1,
               random = ~ gg,
               var.ini = list(gg = 2, resid = 10),
               data = globulus)

Although in most cases the results will not change at all, we encourage to give explicit initial values for variance components. Specially when some estimate can be artifact. This is also useful for checking sensitivity to initial values.

• In globulus, the family (mum) is nested within the provenance (gg)

• This is a matter of codification:

• Otherwise, in both cases we specify the model in the same way:

random = ~ gg + factor(mum)  # note that mum is numeric

• Furthermore, this approach can handle unbalanced and mixed designs

• Standard R notation:

random = ~ gg * factor(mum)

• Not available yet (feature request?)

• Workaround: build the interaction variable manually

• Example: gg and block are crossed factors

dat <- transform(globulus,
                 interaction = factor(gg:bl))
random = ~ gg + bl + interaction

1. Use remlf90() and the globulus dataset to fit
   • a hierarchical model using mum within gg
   • a factorial model using gg and bl
2. Explore the results with summary()
   • is the family (mum) effect relevant?
   • is there any evidence of interaction between gg and bl?

res.h <- remlf90(fixed = phe_X ~ 1,
                 random = ~ factor(mum) + gg,
                 data = globulus)

# Interaction variable
globulus.f <- transform(globulus,
                        gg_bl = factor(gg:bl))

res.f <- remlf90(fixed = phe_X ~ 1,
                 random = ~ gg + bl + gg_bl,
                 data = globulus.f)

• The family effect is not very important, in terms of explained variance
• However, the model is a bit better with it (AIC, logLik)

summary(res)
summary(res.h)

• Looks like the interaction between block and provenance is negligible
• (apart from the fact that it makes no sense at all, and shuld not have been even considered in the first place)
• compare with the model without interaction

summary(res.f)

## result without interaction
res.f0 <- remlf90(fixed  = phe_X ~ 1,
                  random = ~ gg + bl,
                  data = globulus)
paste('AIC:', round(extractAIC(res.f0)),
      'logLik:', round(logLik(res.f0)))

• Random effect at individual level

• Based on a pedigree

• BLUP of Breeding Values from own and relatives' phenotypes

• Represents the additive component of the genetic value

• More general:
  • family effect is a particular case
  • accounts for more than one generation
  • mixed relationships

• More flexible: allows to select individuals within families

res.animal <- remlf90(fixed  = phe_X ~ 1,
                      random = ~ gg,
                      genetic = list(model = 'add_animal', 
                                     pedigree = globulus[, 1:3],
                                     id = 'self'), 
                      data = globulus)

• The pedigree needs to meet certain conditions

• If it does not, breedR automatically completes, recodes and sorts

• If recoding is necessary, breedR issues a warning because you need to be careful when retrieving results

• See this guide for more details

• The residuals of any LMM must be noise

• However, most times there are environmental factors that affect the response

• This causes that observations that are close to each other tend to be more similar that observations that are far away

• This is called spatial autocorrelation

• It may affect both the estimations and their accuracy

• This is why experiments are randomized into spatial blocks

• Include an explicit spatial effect in the model

• I.e., a random effect with a specific covariance structure that reflects the spatial relationship between individuals

• The block effect, is a very particular case:
  • It is designed from the begining, possibly using prior knowledge
  • Introduces independent effects between blocks
  • Most neighbours are within the same block (i.e. share the same effect)

• There seems to remain some intra-block spatial autocorrelation

• nok were (6, 6) by default (see summary())

res.spl99  <- remlf90(fixed  = phe_X ~ 1, random = ~ gg,
                      genetic = gen.globulus,
                      spatial = list(model   = 'splines', 
                                     coord   = globulus[, c('x','y')],
                                     n.knots = c(9, 9)), 
                      data = globulus, method = 'em')

qplot(rho_r, rho_c,
      fill = loglik,
      geom = 'tile',
      data = res.ar1$rho)

res.ar.grid  <- remlf90(fixed  = phe_X ~ gg,
                        genetic = list(model = 'add_animal', 
                                       pedigree = globulus[,1:3],
                                       id = 'self'), 
                        spatial = list(model = 'AR', 
                                       coord = globulus[, c('x','y')],
                                       rho = rho.grid), 
                        data = globulus)
summary(res.ar.grid)

• Optional effect with environmental (rather than genetic) basis

• Modelled as an individual independent random effect that affects neighbouring trees in the same (weighted) way

system.time(
  res.comp <- remlf90(fixed   = phenotype ~ 1,
                      genetic = list(model = 'competition',
                                     pedigree = dat[, 1:3],
                                     id = 'self',
                                     coord = dat[, c('x', 'y')],
                                     competition_decay = 1,
                                     pec = list(present = TRUE)),
                      spatial = list(model = 'AR', 
                                     coord = dat[, c('x', 'y')],
                                     rho   = c(.3, .8)),
                      data = dat,
                      method = 'em')  # AI diverges
  )

##    user  system elapsed 
##  94.458   0.180  94.652

1. Plot the true vs predicted:
   • direct and competition Breeding Values
   • spatial effects
   • pec effects
2. Plot the residuals and their variogram
   • Do you think the residuals are independent?
   • How would you improve the analysis?

## compute the predicted effects for the observations
## by matrix multiplication of the incidence matrix and the BLUPs
pred <- list()
Zd <- model.matrix(res.comp)$'genetic_direct'
pred$direct <- Zd %*% ranef(res.comp)$'genetic_direct'

## Watch out! for the competition effects you need to use the incidence
## matrix of the direct genetic effect, to get their own value.
## Otherwise, you get the predicted effect of the neighbours on each individual.
pred$comp <- Zd %*% ranef(res.comp)$'genetic_competition'
pred$pec  <- as.vector(model.matrix(res.comp)$pec %*% ranef(res.comp)$pec)

comp.pred <-
  rbind(data.frame(Component = 'direct BV', True = dat$BV1, Predicted = pred$direct),
        data.frame(Component = 'competition BV', True = dat$BV2, Predicted = pred$comp),
        data.frame(Component = 'pec', True = dat$pec, Predicted = pred$pec))

ggplot(comp.pred,
       aes(True, Predicted)) +
  geom_point() + 
  geom_abline(int = 0, sl = 1,
              col = 'darkgray') +
  facet_grid(~ Component)

The predition of the Permanent Environmental Competition effect is not precisely great…

plot(res.comp, type = 'resid')
variogram(res.comp)

This additional component allows to introduce a random effect \(\psi\) with arbitrary incidence and covariance matrices \(Z\) and \(\Sigma\):

\[ \begin{aligned} y = & \mu + X \beta + Z \psi + \varepsilon \\ \psi \sim & N(0, \sigma_\psi^2 \Sigma_\psi) \\ \varepsilon \sim & N(0, \sigma_\varepsilon^2) \end{aligned} \]

## Fit a blocks effect using generic
inc.mat <- model.matrix(~ 0 + bl, globulus)
cov.mat <- diag(nlevels(globulus$bl))
res.blg <- remlf90(fixed  = phe_X ~ gg,
                   generic = list(block = list(inc.mat,
                                               cov.mat)),
                   data   = globulus)

## Linear Mixed Model with pedigree and spatial effects fit by AI-REMLF90 ver. 1.110 
##    Data: globulus 
##   AIC  BIC logLik
##  5691 5701  -2844
## 
## Parameters of special components:
## 
## 
## Variance components:
##          Estimated variances   S.E.
## block                  2.592 1.0640
## Residual              15.208 0.6824
## 
## Fixed effects:
##        value   s.e.
## gg.1  13.534 0.6222
## gg.2  14.030 0.8464
## gg.3  16.106 0.5513
## gg.4  11.854 0.6824
## gg.5  15.883 0.5862
## gg.6  10.220 1.3040
## gg.7  13.995 1.0893
## gg.8  15.728 0.5409
## gg.9  16.478 0.5969
## gg.10 12.843 1.1225
## gg.11 16.744 0.6151
## gg.12 17.002 0.8464
## gg.13 16.297 1.0893
## gg.14 14.429 0.4730

• You can predict the Breeding Value of an unmeasured tree

• Or the expected phenotype of a death tree (or an hypothetical scenario)

• Information is gathered from the covariates, the spatial structure and the pedigree

• Simply include the individual in the dataset with the response set as NA

1. Extend the last table to include the predicted value when the phenotype is measured

2. Remove 1/10th of the phenotypes randomly, and predict their expected phenotype
   • Have the parameter estimations changed too much?

3. Compute the Root Mean Square Error (RMSE) of Prediction with respect to the true values

pred.BV.mat <- with(ranef(res.comp), cbind(`genetic_direct`,
                                           `genetic_competition`))

pred.genetic.loo <-
  Zd[rm.idx, ] %*% with(ranef(res.comp),
                        cbind(`genetic_direct`,
                              `genetic_competition`))
valid.pred$Pred.full <- c(Zd[rm.idx, ] %*% pred.BV.mat,
                          fitted(res.comp)[rm.idx])

knitr::kable(round(valid.pred[,c(1,3,2)],
                   2))

rm.idx <- sample(nrow(dat),
                 nrow(dat)/10)
dat.cv <- dat
dat.cv[rm.idx, 'phenotype'] <- NA
## Re-fit the model and build table

true.exp.cv <- with(dat[rm.idx, ], phenotype - resid)
round(sqrt(mean((fitted(res.comp.cv)[rm.idx] - true.exp.cv)^2)), 2)

## [1] 1.5

• We have used simulated data from the metagene software

• If you simulate data, import the results with read.metagene()

• Use several common methods with a metagene object:
  • summary(), plot(), as.data.frame()
• Plus some more specific metagene functions:
  • b.values(), get.ntraits(), ngenerations(), nindividuals(), get.pedigree()
• And specific functions about spatial arrangement:
  • coordinates() extract coordinates
  • sim.spatial() simulates some spatial autocorrelation

• The function breedR.sample.phenotype() simulates datasets from all the model structures available in breedR

• Limitation: only one generation, with random matings of founders

• See ?simulation for details

If you have access to a Linux server through SSH, you can perform computations remotely

• Take advantage of more memory or faster processors

• Parallelize jobs

• Free local resources while fitting models

• See ?remote for details