Initial (co)variance specification with included factor random in multi-trait models, and possibility of working with a repeatability model

63 views
Skip to first unread message

jesus valdes hernandez

unread,
Mar 19, 2020, 6:42:07 AM3/19/20
to breedR

Hello breedR group,

 

I would like you to help me with the following questions and some possible solution for my problem:

 

1. About “Initial (co)variance specification” in multitrait models (pag. 34 of Overview manual):

I try to adapt your example for my case:

 

library(breedR)

 

dat <-

  droplevels(

    douglas[douglas$site == "s3",

            names(douglas)[!grepl("H0[^4]|AN|BR|site", names(douglas))]]

  )

 

initial_covs <- list(

mum = matrix(c(1, .5, .5, 1), nrow = 2),

genetic = 1e3*matrix(c(1, .5, .5, 1), nrow = 2),

residual = diag(2) # no residual covariances

)

 

system.time(

  res <-

    remlf90(

      fixed = cbind(H04, C13) ~ orig,

      random = ~ mum,

      genetic = list(

        model = 'add_animal',

        pedigree = dat[, 1:3],

        id = 'self',

        var.ini = initial_covs$mum,

        var.ini = initial_covs$genetic),

      data = dat,

      var.ini = list(residual = initial_covs$residual)

    )

)

 

But I have the following error:

Error in check_genetic(model = "add_animal", pedigree = list(self = c(8306,  :

  el argumento formal "var.ini" concuerda con múltiples argumentos especificados

Timing stopped at: 0 0 0

 

2.Regarding the specification of the model 'add_animal', the label is implemented to work with a repeatability model?

 

Thank you for your time, and your comments are welcome.

Jesús


jesus valdes hernandez

unread,
Mar 22, 2020, 11:39:34 AM3/22/20
to breedR
Dear Facunado,

I correctly received your email response.
Indeed, I was able to fix the error. Now the correct way would be as follows:

library(breedR)

 

dat <-

  droplevels(

    douglas[douglas$site == "s3",

            names(douglas)[!grepl("H0[^4]|AN|BR|site", names(douglas))]]

  )

 

initial_covs <- list(

mum = 1e3*matrix(c(1, .5, .5, 1), nrow = 2),

genetic = 1e3*matrix(c(1, .5, .5, 1), nrow = 2),

residual = diag(2) # no residual covariances

)

 

system.time(

  res <-

    remlf90(

      fixed = cbind(H04, C13) ~ orig,

      random = ~ mum,

      genetic = list(

        model = 'add_animal',

        pedigree = dat[, 1:3],

        id = 'self',

        var.ini = initial_covs$genetic),

      data = dat,

      var.ini = list(mum = initial_covs$mum, residual = initial_covs$residual)

    )

)


On the other hand, I only have 2 questions regarding:
How can you declare more than one random effect in the breedR syntax (e.g., animal and permanent environmental effect)?
Just as I would declare the argument var.ini considering this particular case (i.e. 2 random effects, the genetic part and the residual part).

Thanks for your time,
Best regards,
Jesús

Facundo Muñoz

unread,
Mar 22, 2020, 2:04:59 PM3/22/20
to bre...@googlegroups.com

Dear Jesús,

How is is your model for the permanent environmental effect? Isn't it covered by one of the spatial effects (e.g. the autoregressive model)?

In that case, the initial variance can be specified within the `spatial` component (just like the genetic). See examples in the Overview tutorial, or in the help for remlf90 (`?remlf90`).

ƒacu.-

--
Report issues at https://github.com/famuvie/breedR/issues
---
You received this message because you are subscribed to the Google Groups "breedR" group.
To unsubscribe from this group and stop receiving emails from it, send an email to breedr+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/breedr/461f7b43-7441-416e-a4b6-5d164a794c6d%40googlegroups.com.

jesus valdes hernandez

unread,
Mar 27, 2020, 12:58:38 PM3/27/20
to breedR

Hello dear ƒacu.-


I have prepared a simulated example to explain what I want to learn and thus obtain your recommendations and suggestions.


1.How is is your model for the permanent environmental effect?

A1:

#Specifying a pedigree (example with 10 individuals and repated measures):

self <- as.numeric(paste("1",sep="0", c(rep("1",5),rep("2",1),rep("3",3),rep("2",2),rep("4",1),rep("5",1),rep("3",1),rep("4",1),rep("5",3),rep("6",1),rep("7",2),rep("8",1),rep("9",2),rep("8",1),rep("9",1),rep("6",2),rep("10",1),rep("8",1),rep("4",1),rep("10",1),rep("11",4))))

dad <- as.numeric(c(0,0,0,0,0,paste("0", sep="0",c(2,1,1,1,2,2,1,2,1,1,2,2,2,4,4,4,6,6,6,6,6,4,4,7,6,1,7,7,7,7,7))))

mum <- as.numeric(c(0,0,0,0,0,paste("0",sep="0",c(3,3,3,3,3,3,5,5,3,5,5,5,5,8,9,9,9,8,8,9,8,8,8,9,9,5,9,10,10,10,10))))

 

#Factors and phenotypes:

sex <- as.factor(c(rep("1",5),rep("2",1),rep("1",3),rep("2",2),rep("2",1),rep("1",1),rep("1",1),rep("2",1),rep("1",3),rep("2",1),rep("1",2),rep("2",1),rep("1",2),rep("2",1),rep("1",1),rep("2",2),rep("1",1),rep("2",1),rep("2",1),rep("1",1),rep("2",4)))

phen1 = seq(min(5.0),max(31.5),length=36)

phen2 = seq(min(1.6),max(3.85),length=36)

animal <- as.factor(c(rep("1",5),rep("2",1),rep("3",3),rep("2",2),rep("4",1),rep("5",1),rep("3",1),rep("4",1),rep("5",3),rep("6",1),rep("7",2),rep("8",1),rep("9",2),rep("8",1),rep("9",1),rep("6",2),rep("10",1),rep("8",1),rep("4",1),rep("10",1),rep("11",4)))

pef = as.data.frame(rank(animal))

pef[] <- lapply(pef, as.integer)

pef <- pef$`rank(animal)`

 

data = cbind.data.frame(self,dad,mum,sex,animal,x,y,pef,phen1,phen2)


library(breedR)

#Fitting an animal model (Multiple traits and two random effects):

#Initial (co)variance specification:

initial_covs <- list(

  genetic = 1e3*matrix(c(1, .5, .5, 1), nrow = 2),

  animal = 1e3*matrix(c(7.544e-02,6.391e-03,6.391e-03,9.608e-04), nrow = 2),

  pef = 1e3*matrix(c(7.544e-02,6.391e-03,6.391e-03,9.608e-04), nrow = 2),

  residual = diag(2)

)

# pef is the permanent environment effect

res.animal <- remlf90(fixed = cbind(phen1,phen2) ~ sex,

                      random = ~ animal + pef,

                      genetic = list(model = 'add_animal',

                                     pedigree = data[, 1:3],

                                     id = 'self',

                                     var.ini = initial_covs$genetic),

                      data = data,

                      var.ini = list(animal = initial_covs$animal, pef=initial_covs$pef, residual = initial_covs$residual))

summary(res.animal)


2.Isn't it covered by one of the spatial effects (e.g. the autoregressive model)?

A2:

#The genetic component:

gen.part <- list(model = 'add_animal',

                 pedigree = data[, 1:3],

                 id = 'self',

                 method = 'em')

#2.1.The blocks model:Animal-spatial model: results

res.blk <- remlf90(fixed = cbind(phen1,phen2) ~ sex,

                   random = ~ animal,

                   genetic = gen.part,

                   spatial = list(model = 'blocks',

                                  coord = data[, c('x', 'y')],

                                  id = 'pef'),

                   data = data)

summary(res.blk)

#2.2.Autoregressive model:

res.ar1 <- remlf90(fixed = cbind(phen1,phen2) ~ sex,

                      random = ~ animal,

                      genetic = gen.part,

                   spatial = list(model = 'AR',

                                  coord = data[, c('x','y')],

                                  rho   = c(.85, .8)),

                                  data = data)

summary(res.ar1)

 

Comments and questions:

Based on help for ?breedR.options and ?remlf90


There are several issues that I was not able to understand well, especially for examples 2.1 and 2.2. In fact:

a)      a) I have not been able to find how it calculates the “x” and “y” coordinates. 

       -How would you calculate them in my case?


b)      b)Also, in the help I am not clear how to declare in the 2.1 syntaxis. and 2.1 the var.ini argument. This is because I have two random factors (animal and pef). Indeed, it is very important to be able to declare "Initial (co) variance specification" because I have noticed that in this way the estimates are correct, and for example it does not come out “NA in AIC BIC and logLik criteria”.

S     -So, how would the var.ini argument be declared for these specific cases?

 

Finally, sorry for the length of the post.


I look forward to your response and help.

Best regards,

Jesús

To unsubscribe from this group and stop receiving emails from it, send an email to bre...@googlegroups.com.

Facundo Muñoz

unread,
Apr 5, 2020, 6:48:22 AM4/5/20
to bre...@googlegroups.com

Dear Jesús,

Thank you for your patience: I'm not having much time these days (less than usual!). I know you understand.

Here are somme comments on your example. I hope they help. Don't hesitate to continue the discussion if you have further questions.


1. x and y were not defined when constructing `data`. Just to be able to
run your code I used:

x <- rep(1:6, each = 6)
y <- rep(1:6, times = 6)


2. In the first model (res.animal), animal and pef are not identifiable, since they both have
one estimated value per animal, by construction of pef. Even if their values
in the dataset are different, there is a one to one correspondence between them.
unique(data[, c("animal", "pef")])

Consider the data values of random effects as mere "labels".

Besides, it does not make much sense to associate and "environmental" effect
to the individuals. Suppose that the first five individuals are in different
locations. Why would they share the same pef? Unless they are located in the
same plot, in which case it is impossible to separate both effects.

Moreover, these animal and pef effects are also somewhat confounded with the
additive genetic effect, since again there is a one to one correspondence
between their values/labels:
 
identical(as.numeric(levels(animal)[animal]), as.numeric(factor(self)))  # TRUE

The difference is that the genetic effects are not independent from
one another (they are related through the pedigree). This allows to
separate and quantify to some extent the contribution of the genetic effect
to that of the `animal` or the `pef`, but not both!!.


3. In the subsequent models, `method = 'em'` should not be part of the
genetic term (`gen.part`). It should go outside, as an argument to
`remlf90()` function. Note that you are using AI-REML in all your models.

4. The model `res.blk` is identical (mathematically) to `res.animal`,
in the sense that the spatial blocks model (from `res.blk`) is equivalent to specifying
the block variable (`pef`) in the `random` component (as in `res.animal`).
So, the same comments apply.

In practice, some negligible differences in the results might be due to
not specifying initial covariances (as in `res.animal`).


5. In contrast, the model `res.ar1` fits an spatial model based on the
locations (coordinates x and y) of the individuals rather than on `pef`
(which is associated to the genotypes, not the locations).
In this sence, it is more properly an environmental effect since it
will capture variations due to the spatial location.


Stay safe. Best wishes.

ƒacu.-

To unsubscribe from this group and stop receiving emails from it, send an email to breedr+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/breedr/3e7b2173-5910-44f4-ba26-00385ffd0764%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages