4.6 Dominance

Here we can make use of the dominance relatedness matrices that can be generated in the nadiv package

data(warcolak)

pedD <- warcolak[,1:3]

Dmats <- makeD(pedD)

## starting to make D....done 
## starting to invert D....done

Dmat <- Dmats$D

We can input this into simulate_population() with the cov_str argument:

ds <- data.frame(animal = pedD[,1], animalD = pedD[,1])

squid_data <- simulate_population(
  data_structure = ds, 
 
  pedigree = list(animal=pedD[,1:3]),
  cov_str = list(animalD=Dmat),

  parameters =list(
    animal = list(
      vcov = 0.2
    ),
    animalD = list(
      vcov = 0.2
    ),
    residual = list(
      vcov = 0.5
    )
  )
)

data <- get_population_data(squid_data)

We can then run a model to check the simulation, using {gremlin}, with inverse A and D matrices:

Dinv <- Dmats$Dinv
Ainv_dom <- makeAinv(pedD)$Ainv

mod_dom <- gremlin(y~1, random=~ animal + animalD,data=data,ginverse=list(animal=Ainv_dom,animalD=Dinv))

## gremlin started:      11:24:45 
##   1 of max 20        lL:-5762.227161     took 0.0061 sec.
##   2 of max 20        lL:-2494.762867     took 0.0049 sec.
##   3 of max 20        lL:-2312.124280     took 0.0049 sec.
##   4 of max 20        lL:-2299.455179     took 0.0048 sec.
##   5 of max 20        lL:-2299.355390     took 0.0048 sec.
##   6 of max 20        lL:-2299.355382     took 0.0048 sec.
##   7 of max 20        lL:-2299.355382     took 0.0048 sec.
## 
## ***  REML converged  ***
## 
## gremlin ended:        11:24:45

summary(mod_dom)

## 
##  Linear mixed model fit by REML [' gremlin ']
##  REML log-likelihood: -2299.355 
##   lambda: FALSE 
## 
##  elapsed time for model: 0.0863 
## 
##  Scaled residuals:
##      Min       1Q   Median       3Q      Max 
## -2.92188 -0.54659  0.02061  0.55849  3.09950 
## 
##  (co)variance parameters: ~animal + animalD 
##          rcov: ~units 
##           Estimate Std. Error
## G.animal    0.2068    0.02422
## G.animalD   0.1588    0.04980
## ResVar1     0.5437    0.04572
## 
##  (co)variance parameter sampling correlations:
##           G.animal G.animalD  ResVar1
## G.animal   1.00000   -0.2807 -0.01946
## G.animalD -0.28069    1.0000 -0.89553
## ResVar1   -0.01946   -0.8955  1.00000
## 
##  Fixed effects: y ~ 1 
##             Estimate Std. Error z value
## (Intercept) -0.02196    0.02418 -0.9084