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
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)
)
)
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.
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.
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.