Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Different results with different classes of the 'id' variable #76

Open
nalnahhas opened this issue Dec 21, 2016 · 4 comments
Open

Different results with different classes of the 'id' variable #76

nalnahhas opened this issue Dec 21, 2016 · 4 comments

Comments

@nalnahhas
Copy link

Hi there,
I recently tried to fit an animal model using the remlf90() function. My model was simple and contained 4 fixed effects, 1 random non-genetic effect and the genetic additive effect (pedigree). I compared the results (h2 + se) to those of BLUPF90 (airemlf90) and they were the same as they should be. Then, I changed the class of the 'id' variable in the genetic part of the model from integer to factor and I re-ran the model. The h2 was considerabley different from that I got when the 'id' was of class integer.

dat399$animal <- as.interger(dat399$animal) # h2 = 0.44, se = 0.012 (correct)
dat399$animal <- factor(dat399$animal) # h2 = 0.08, se = 0.006 (wrong)
dat399$animal <- as.interger(as.character(dat399$animal)) # h2 = 0.44, se = 0.012 (correct again)

So, is this normal, should the 'id' part of the genetic effect be always coded as integer or there is a bug that needs to be corrected?

Many thanks,

Nabeel.

@famuvie
Copy link
Owner

famuvie commented Dec 21, 2016

Hi Nabeel.
Yes. The variables encoding individuals (i.e., id and progenitors) should be integers.
However, the pedigree-building function should have raisen an error whenever the user tries with a factor.
How did you specify the pedigree in the model?

Thanks for your report and help.

@nalnahhas
Copy link
Author

Hi Facundo,
I specified the model using coventional R code:

mod.anim.SHI <- remlf90(fixed = SHI ~ lot + date + plato + pos,
random = ~ repet,
genetic = list(model = 'add_animal',
pedigree = ped399[, c(1:3)],
id = 'animal'),
data = dat399)

Many thanks.

@nalnahhas
Copy link
Author

Facundo,
On a different matter, when you read text files using the 'readr' package, the R object containing the data will have both classes, tibble and data.frame. When using this kind of object with the Breedr package, an error message is sent to the console and you are forced to read your text files using functions from the base in order to get a data.frame object that you can use with Breedr without any problem. So, the breedr package is not compatible with tibbles or I'm doing something wrong while reading my text files into R?

Best regards,

Nabeel.

@famuvie
Copy link
Owner

famuvie commented Jan 9, 2017

Dear Nabeel,
thanks for your feedback.

  1. Regarding the issue of the variable type of the animal id, note that in the genetic component, the pedigree is taken from ped399, where the variable animal is presumably integer or numeric. However, if you change the corresponding variable in dat399 as you have been doing, this breaks the correspondence between the animal codes in the pedigree and the dataset.

  2. Concerning tibbles, I had not look at it so far, but it seems to work for me:

library(tibble)
summary(
  remlf90(y ~ 1, data = tibble(y = rnorm(10)))
)
#> No specification of initial variances.
#>      Using default value of 1 for all variance components.
#>      See ?breedR.getOption.
#> Linear Mixed Model with pedigree and spatial effects fit by AI-REMLF90 ver. 1.122 
#>    Data: tibble(y = rnorm(10)) 
#>    AIC   BIC logLik
#>  29.83 30.13 -13.92
#> 
#> Parameters of special components:
#> 
#> 
#> Variance components:
#>          Estimated variances   S.E.
#> Residual                1.06 0.4999
#> 
#> Fixed effects:
#>             value   s.e.
#> Intercept 0.15672 0.3257

Which error message did you get exactly?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants