Using Bootstrap estimation for a Spline model
Pamela Villalba
Hello to all breedR users,
I'm trying to use a bootstrap estimation for a spline model. Because I want to estimate the s.e. for the model variances (given that the "em" method can't estimate them). Has anyone already tried to do this?
I’m using the globulus data from the package as an example, but I’m having problems with it.
Would someone be so kind to give me a hand with this?
Thanx,
Pamela
Facundo Muñoz
Dear Pamela,
have you checked the guide on bootstrap computation of heritability [1]?
You can easily adapt the function `sim_target()` in step 3 to recover the model variances (instead of the heritability as it is done in the guide).
I hope it helps.
ƒacu.-
[1]
https://github.com/famuvie/breedR/wiki/Heritability#case-3-using-bootstrap-estimation
--
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 post to this group, send email to bre...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/breedr/99014698-e733-4db8-90e8-ba8d3e5c3e45%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Pamela Villalba
Facundo Muñoz
Indeed, the sampling function is always `breedR.sample.phenotype()`.
Even if you are simulating a splines model, the splines are only a component of the model.
Also, note that in function sim_target() you need to fit the model using `data = dat` to fit the model
to the simulated data. Otherwise you will get always the same results.
Best
ƒacu.-
----- Mail original -----
De: "Pamela Villalba" <pvill...@gmail.com>
À: "breedR" <bre...@googlegroups.com>
Envoyé: Mardi 12 Février 2019 17:01:09
Objet: Re: [breedR] Using Bootstrap estimation for a Spline model
> To post to this group, send email to bre...@googlegroups.com <javascript:>
> .
> https://groups.google.com/d/msgid/breedr/99014698-e733-4db8-90e8-ba8d3e5c3e45%40googlegroups.com
> .
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 post to this group, send email to bre...@googlegroups.com.
Pamela Villalba
Hello Facundo,
Is it possible to simulate fixed effect values such as genetic groups (gg in globulus dataset) with the "breedR.sample.phenotype ()" function of breedR?
Because I made the corresponding changes that you mentioned previously, but now I have found a new error and I do not know if it is possible to solve it. The error is within the fixed effects values in the "dat" object. I show you the script and the error that I have next:
|
> dat <- breedR.sample.phenotype( + fixed = phe_X ~ gg, + random = list(fam = list(nlevels = 63, + sigma2 = res$var['fam', 1])), + genetic = list(model = 'add_animal', + pedigree = globulus[,1:3], + id = 'self', + var.ini = res$var['genetic',1]), + spatial = list(model = 'splines', + coord = globulus[, c('x','y')], + n.knots = c(7,7), + var.ini = res$var['spatial',1]), + residual.variance = res$var['Residual', 1], + N = nrow(model.frame(res), + method = 'em') + )
Error in X %*% fixed : requires numeric/complex matrix/vector arguments |
|
|
|
Since I want to get the blue values for each gg; I tried to include the fixed effect value (for the globulus dataset: 14 genetic groups) in the function as shown below:
> X=model.matrix(~-1+globulus$gg)
> dat <- breedR.sample.phenotype(
+ fixed = X,
+ random = list(bl = list(nlevels = 15,
+ sigma2 = res$var['spatial', 1]),
+ fam = list(nlevels = 63,
+ sigma2 = res$var['genetic', 1])),
+ residual.variance = res$var['Residual', 1],
+ N = nrow(model.frame(res))
+ )
Error in X %*% fixed : non-conformable arguments
However, I still find an error in the fixed parameter.
So it’s possible to simulate this type of fixed values with the "breedR.sample.phenotype ()" function? Especially if it isn’t a regression coefficient?
Thank you!!
Facundo Muñoz
Dear Pamela,
sorry for the late reply.
It is not possible to
do it within the breedR.sample.phenotype()
function, but it is very easy to work it around and do it
manually. Indeed, this function is not fully complete, but is
focused on the more complex terms.
The fixed argument in this function requires a (possibly named) numeric vector. Not a formula (as in your first example), nor a matrix (as in your second example). The numbers in this vectors are interpreted as regression coefficients of random uniform covariates. This is not what you want.
Furthermore, the
syntax in your first example is not accurate: you wrote the genetic
and sptaial
arguments as they are for the fitting function remlf90(),
not as they are for the simulator breedR.sample.phenotype().
See the help page.
So, the work around is :
-
Simulate a dataset with only a global intercept in the fixed part and all the random components you need
-
Manually add columns with categorical (factor) covariates (like genetic groups) to the dataset
-
Manually modify the phenotype according to values of the categories
Here is an example:
library(breedR)
#> Loading required package: sp
## 1. Simulate dataset
dat <- breedR.sample.phenotype(
fixed = c(mu = 1),
genetic = list(
model = 'add_animal',
Nparents = c(10, 10),
sigma2_a = 1,
check.factorial = FALSE
),
spatial = list(
model = 'splines',
grid.size = c(10, 10),
n.knots = c(7, 7),
sigma2_s = 1
),
residual.variance = 1
)
## Remove founders
dat <- dat[!is.na(dat$sire),]
## 2. Add categorical covariate
## assign categories randomly
dat$gg <- sample(factor(letters[1:14]), nrow(dat), replace = TRUE)
## 3. Modify phenotype according to categories
## vector of deviations from global mean for each category
gg_values <- setNames(rnorm(nlevels(dat$gg)), levels(dat$gg))
dat$phenotype <- dat$phenotype + gg_values[dat$gg]
I hope it helps
ƒacu.-
To view this discussion on the web visit https://groups.google.com/d/msgid/breedr/045afe4c-c9de-4289-9a30-5b09353aabc6%40googlegroups.com.