Diagnosis of spatial and competiton models

10-11 Dec. 2015, breedR version: 0.11

Douglas datset

Description of Douglas datset

focus on Site 3

data("douglas")
str(douglas)

'data.frame':   9630 obs. of  15 variables:
 $ self : num  135 136 137 138 139 140 141 142 143 144 ...
 $ dad  : num  41 41 41 41 41 41 41 41 41 41 ...
 $ mum  : num  21 21 21 21 21 21 21 21 21 21 ...
 $ orig : Factor w/ 11 levels "pA","pB","pC",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ site : Factor w/ 3 levels "s1","s2","s3": 1 1 1 1 1 1 1 1 1 1 ...
 $ block: Factor w/ 127 levels "s1:1","s1:2",..: 11 44 24 28 13 35 8 40 3 15 ...
 $ x    : num  6 27 45 57 57 60 63 66 66 75 ...
 $ y    : num  81 135 90 45 327 474 450 21 234 483 ...
 $ H02  : int  NA NA NA NA NA NA NA NA NA NA ...
 $ H03  : int  NA NA NA NA NA NA NA NA NA NA ...
 $ H04  : int  NA NA NA NA NA NA NA NA NA NA ...
 $ H05  : int  634 581 611 370 721 488 574 498 528 620 ...
 $ C13  : int  586 474 715 372 665 558 490 372 527 612 ...
 $ AN   : Factor w/ 5 levels "1","2","3","4",..: NA NA NA NA NA NA NA NA NA NA ...
 $ BR   : Factor w/ 5 levels "1","2","3","4",..: NA NA NA NA NA NA NA NA NA NA ...

dat <- droplevels(subset(douglas, site == "s3"))

Circumference 2013

ggplot(dat, aes(C13)) + geom_histogram()

ggplot(dat, aes(x, y, fill = C13)) +
  geom_raster(show_guide = FALSE) +
  coord_fixed() + 
  scale_fill_viridis()

Variogram of phenotype

Slight autocorrelation, more apparent across rows

variogram(coord = dat[, c('x','y')], z = dat$C13)

Basic model

family + block design effects

Since there is no explicit genetic effect, we give the formula for the heritability

dat$fam <- factor(dat$mum)
h2.fml <- "4*G_3_3_1_1/(G_2_2_1_1+G_3_3_1_1+R_1_1)"

res.base <- remlf90(
  fixed = C13 ~ orig,
  random = ~ block + fam,
  progsf90.options = paste("se_covar_function h2", h2.fml),
  dat = dat)

Linear Mixed Model with pedigree and spatial effects fit by AI-REMLF90 ver. 1.122 
   Data: dat 
   AIC   BIC logLik
 17987 18003  -8991

Parameters of special components:


Variance components:
         Estimated variances S.E.
block                   1087  387
fam                     1774  604
Residual               19914  774

   Estimate   S.E.
h2   0.3095 0.1004

Fixed effects:
         value   s.e.
orig.pA 471.86 12.470
orig.pB 501.58 19.925
orig.pC 435.96 27.087
orig.pF 444.08 14.372
orig.pG 378.28 50.751
orig.pH 389.80 46.890
orig.pI 409.80 47.066
orig.pJ 416.98 46.660
orig.pK 445.68 46.509

Individual-tree model

additive-genetic effect + block design effects

genetic component -> heritability for free
use a proper blocks spatial model (convenient for plotting)

res.blk <- remlf90(
  fixed = C13 ~ orig,
  genetic = list(model = 'add_animal',
                 pedigree = dat[, c('self','dad','mum')],
                 id = 'self'),
  spatial = list(model = 'blocks',
                 coord = dat[, c('x','y')],
                 id = "block"),
  dat = dat)

Linear Mixed Model with pedigree and spatial effects fit by AI-REMLF90 ver. 1.122 
   Data: dat 
   AIC   BIC logLik
 17984 18000  -8989

Parameters of special components:
spatial: n.blocks: 34

Variance components:
         Estimated variances S.E.
genetic                 5477 1727
spatial                 1107  391
Residual               16235 1384

             Estimate    S.E.
Heritability   0.2389 0.07069

Fixed effects:
         value   s.e.
orig.pA 466.96 14.430
orig.pB 493.67 20.055
orig.pC 435.60 25.267
orig.pF 444.01 13.464
orig.pG 378.21 46.773
orig.pH 389.83 42.457
orig.pI 409.72 42.655
orig.pJ 417.02 42.198
orig.pK 445.65 42.027

Splines model

additive-genetic effect + bidimensional Bsplines

Since we need to use method = em with the splines, we loose the automatic computation of heritability,
Excercise: implement the bootstrap strategy explained in the vignette('Heritability')

res.spl <- remlf90(
  fixed = C13 ~ orig,
  genetic = list(model = 'add_animal',
                 pedigree = dat[, c('self','dad','mum')],
                 id = 'self'),
  spatial = list(model = 'splines',
                 coord = dat[, c('x','y')]),
  dat = dat,
  method = 'em')

Linear Mixed Model with pedigree and spatial effects fit by REMLF90 ver. 1.79 
   Data: dat 
   AIC     BIC logLik
 17963 unknown  -8978

Parameters of special components:
spatial: n.knots: 13 13

Variance components:
         Estimated variances
genetic                 5430
spatial                 1167
Residual               16210

Fixed effects:
         value   s.e.
orig.pA 471.81 19.360
orig.pB 496.26 23.797
orig.pC 441.95 28.341
orig.pF 448.04 18.618
orig.pG 384.14 48.376
orig.pH 394.08 44.184
orig.pI 412.29 44.410
orig.pJ 421.61 43.963
orig.pK 448.73 43.842

Autoregressive model

additive-genetic effect + AR1xAR1

Excercises:

compare with the model including the blocks effects also
find optimal values for rho

res.ar1 <- remlf90(
  fixed = C13 ~ orig,
  # random = ~ block,
  genetic = list(model = 'add_animal',
                 pedigree = dat[, c('self','dad','mum')],
                 id = 'self'),
  spatial = list(model = 'AR',
                 coord = dat[, c('x','y')],
                 rho = c(.8,.8)),
  dat = dat)

Linear Mixed Model with pedigree and spatial effects fit by AI-REMLF90 ver. 1.122 
   Data: dat 
   AIC   BIC logLik
 17981 17997  -8987

Parameters of special components:
spatial: rho: 0.8 0.8

Variance components:
         Estimated variances   S.E.
genetic                 5613 1745.0
spatial                 1342  438.4
Residual               15611 1409.3

             Estimate    S.E.
Heritability   0.2475 0.07187

Fixed effects:
         value   s.e.
orig.pA 467.27 14.735
orig.pB 493.71 20.310
orig.pC 437.33 25.581
orig.pF 443.79 13.720
orig.pG 381.11 47.239
orig.pH 390.88 42.919
orig.pI 408.63 43.114
orig.pJ 417.51 42.669
orig.pK 445.64 42.527

Comparison of residuals from different models

coordinates(res.base) <- dat[, c('x', 'y')]
compare.plots(    # preserve scale
  list(`Family/blocks`      = plot(res.base, 'residuals'),
       `Individual/blocks`  = plot(res.blk, 'residuals'),
       `Individual/splines` = plot(res.spl, 'residuals'),
       `Individual/AR1`     = plot(res.ar1, 'residuals'))
)

Comparison of spatial components

compare.plots(
  list(Blocks = plot(res.blk, type = 'spatial'),
       Splines = plot(res.spl, type = 'spatial'),
       AR1xAR1 = plot(res.ar1, type = 'spatial'))
)

Competition models

Competition genetic model + block spatial model

Excercise: compare with using a decay \(\alpha = 1\)

res.comp <- remlf90(
  fixed = C13 ~ orig,
  genetic = list(
    model = c('comp'),
    pedigree = dat[, c('self','dad','mum')],
    id = 'self',
    coord = dat[, c('x', 'y')],
    competition_decay = 2, # IC decay 1/distance
    pec = TRUE), # envirmonetal compettion effect
  spatial = list(
    model = 'blocks',
    coord = dat[, c('x','y')],
    id = "block"),
  data = dat,
  method = 'em'
)

summary(res.comp)

Linear Mixed Model with pedigree and spatial effects fit by REMLF90 ver. 1.79 
   Data: dat 
   AIC     BIC logLik
 17984 unknown  -8986

Parameters of special components:
spatial: n.blocks: 34

Variance components:
$genetic
                    genetic_direct genetic_competition
genetic_direct                5481             -1240.0
genetic_competition          -1240               400.7

$pec
      [,1]
[1,] 401.4

$spatial
     [,1]
[1,] 1111

$Residual
      [,1]
[1,] 15450


Fixed effects:
         value   s.e.
orig.pA 465.53 12.287
orig.pB 487.29 19.126
orig.pC 437.31 24.652
orig.pF 448.89 12.970
orig.pG 360.13 45.398
orig.pH 382.20 39.633
orig.pI 424.90 40.445
orig.pJ 425.57 39.400
orig.pK 428.71 39.303

direct-competition genetic correlation

(S <- res.comp$var$genetic)

                    genetic_direct genetic_competition
genetic_direct                5481             -1240.0
genetic_competition          -1240               400.7

SD <- sqrt(diag(1/diag(S)))
(SD %*% S %*% SD)[1,2]

[1] -0.8367238

Competition genetic model + AR1xAR1 spatial model

res.comp.ar1 <- remlf90(
  fixed = C13 ~ orig,
  genetic = list(
    model = c('comp'),
    pedigree = dat[, c('self','dad','mum')],
    id = 'self',
    coord = dat[, c('x', 'y')],
    competition_decay = 2, # IC decay 1/distance
    pec = list(present = TRUE)), #envirmonetal compettion effect
  spatial = list(
    model = 'AR',
    coord = dat[, c('x','y')],
    rho = c(.8,.8)),
  data = dat,
  method = 'em'
)

summary(res.comp.ar1)

Linear Mixed Model with pedigree and spatial effects fit by REMLF90 ver. 1.79 
   Data: dat 
   AIC     BIC logLik
 17978 unknown  -8983

Parameters of special components:
spatial: rho: 0.8 0.8

Variance components:
$genetic
                    genetic_direct genetic_competition
genetic_direct                6058             -1553.0
genetic_competition          -1553               504.5

$pec
      [,1]
[1,] 211.3

$spatial
     [,1]
[1,] 1500

$Residual
      [,1]
[1,] 14440


Fixed effects:
         value   s.e.
orig.pA 465.57 12.646
orig.pB 485.89 19.570
orig.pC 439.73 25.429
orig.pF 449.19 13.438
orig.pG 358.90 46.738
orig.pH 382.45 40.506
orig.pI 427.66 41.408
orig.pJ 429.96 40.269
orig.pK 423.57 40.280

Competition genetic model + splines spatial model

res.comp.spl <- remlf90(
  fixed = C13 ~ orig,
  genetic = list(
    model = c('comp'),
    pedigree = dat[, c('self','dad','mum')],
    id = 'self',
    coord = dat[, c('x', 'y')],
    competition_decay = 2, # IC decay 1/distance
    pec = list(present = TRUE)), #envirmonetal compettion effect
  spatial = list(
    model = 'splines',
    coord = dat[, c('x','y')]),
  data = dat,
  method = 'em'
)

summary(res.comp.spl)

Linear Mixed Model with pedigree and spatial effects fit by REMLF90 ver. 1.79 
   Data: dat 
   AIC     BIC logLik
 17962 unknown  -8975

Parameters of special components:
spatial: n.knots: 13 13

Variance components:
$genetic
                    genetic_direct genetic_competition
genetic_direct                5595             -1296.0
genetic_competition          -1296               449.3

$pec
     [,1]
[1,]  423

$spatial
     [,1]
[1,] 1163

$Residual
      [,1]
[1,] 15250


Fixed effects:
         value   s.e.
orig.pA 470.92 17.915
orig.pB 489.03 23.082
orig.pC 444.47 27.900
orig.pF 452.44 18.270
orig.pG 367.00 47.382
orig.pH 385.10 41.714
orig.pI 427.34 42.537
orig.pJ 430.68 41.528
orig.pK 428.89 41.531

Basic model comparison

log-likelihood and AIC

ml <- list(`no comp.` = res.blk, `comp+blk` = res.comp, 
           `comp+ar1` = res.comp.ar1, `comp+spl` = res.comp.spl)
kable(cbind(loglik = sapply(ml, logLik),
            AIC = sapply(ml, function(x) x$fit$AIC)), digits=2)

	loglik	AIC
no comp.	-8989.05	17984.09
comp+blk	-8986.18	17984.36
comp+ar1	-8983.16	17978.31
comp+spl	-8975.20	17962.39

Excercise

heritability with parallel computing

Use the bootstrap strategy to compute the heritability for both the direct and competition genetic effects.
Use the model with blocks spatial effects (res.comp)
Sample \(N=10\) heritability values each
Then we will join the samples