Initial (co)variance specification with included factor random in multi-trait models, and possibility of working with a repeatability model
jesus valdes hernandez
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
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)
)
)
Facundo Muñoz
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
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
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.