Dear all,
I am using MCMCglmm function to construct bivariate animal models of
bustard sperm production according to age-classes.
The problem is that the models can stochastically crash before the end
of the run (at 2000 iterations or 120000 or other) or can finish. For
the model which does not finish, R returns different errors as:
-Mixed model equations singular: use a (stronger) prior
-ill-conditioned G/R structure: use proper priors if you haven't or
rescale data if you have
For the models which reach the end, the estimations of genetic
additive variance appear quite good (nice bell shaped posterior
disctribution).
The problem still remains when I remove the animal term.
When I run univariate models, it works fine and the posteriors
distributions look very good.
Strangely, the more data I have, the more models crash (the largest
amount of data I have is 65000 data for 2400 individuals for one model).
The model looks like:
priorExp<-list(G=list(G1=list(V=diag(2), nu=2,
alpha.mu=rep(0,2),alpha.V=diag(2)*100000),
G2=list(V=diag(2), nu=2, alpha.mu=rep(0,2),alpha.V=diag(2)*100000),
G3=list(V=diag(2), nu=2, alpha.mu=rep(0,2),alpha.V=diag(2)*100000)),
R=list(V=diag(2), nu=2))
spz<-MCMCglmm(cbind(age1_2,age5_6)~trait-1 + trait:tse+
trait:I(tse^2)+ trait:joe + trait:I(joe^2),
random=~us(trait):animal+us(trait):ID+us(trait):annee ,
rcov=~us(trait):units,nitt=150000, thin=1000, burnin=100000,
prior=priorExp, verbose=TRUE, pedigree=ped,
family=c("gaussian","gaussian"), data=dat)
For the fixed effects : I use 4 continuous parameters as correction
for each trait
For the random effects: I use, individuals, years and animal parameters
I have also tried more informative prior (as described in WAMWIKI) but
the problem was the same.
To give you an example :
Because of computing limitation, I use multi-chain process. I run
several times the same model (as above) and concatenate results (same
prior,same burning, same thin and random seed) to obtain at least 1000
estimates (50 estimates by model). In this context, I ran 50
bivariable models with the age-class age1_2 and the age-class age5_6
but only 9 models of the 50 models reached the end.
When we look fixed parameters estimates (estimate are binded for the
nine models : http://ubuntuone.com/3Gi8GwjcRk3P01MxJp2qLe ), we can
see that the estimates are really close to 0. Could it be the problem?
When we look ramdom parameters estimates (estimate are binded for the
nine models : http://ubuntuone.com/42oaP9euG1m2LNipMawcHX ), the
residual estimates do not look very good. Could it be the problem?
Last thing, if I try to add a cubic effect, all the models crash (same
error than before or memory mapped error).
I really do not know where the problem comes from. Do you have an idea?
Thanks
--
Stephane Chantepie
CNRS Phd candidate
Mus?um national d'Histoire naturelle
55 rue Buffon
75005 paris
E-mail : chantepie_at_mnhn.fr
--
Bivariate animal models with both "ill-conditioned G/R structure" and "Mixed model equations singular" errors
5 messages · Jarrod Hadfield, Stephane Chantepie
Dear Stephane, When you say crash do you mean crash in the sense of a segfault or in the sense that it stops with the errors: -Mixed model equations singular: use a (stronger) prior -ill-conditioned G/R structure: use proper priors if you haven't or rescale data if you have If the latter, it may just require a rescaling of your continuous covariates by using something like scale(). If the former, it would be good for me to have a reproducible example as it means there is a bug. Cheers, Jarrod uoting chantepie at mnhn.fr on Fri, 09 Mar 2012 10:33:42 +0100:
Dear all,
I am using MCMCglmm function to construct bivariate animal models of
bustard sperm production according to age-classes.
The problem is that the models can stochastically crash before the
end of the run (at 2000 iterations or 120000 or other) or can
finish. For the model which does not finish, R returns different
errors as:
-Mixed model equations singular: use a (stronger) prior
-ill-conditioned G/R structure: use proper priors if you haven't or
rescale data if you have
For the models which reach the end, the estimations of genetic
additive variance appear quite good (nice bell shaped posterior
disctribution).
The problem still remains when I remove the animal term.
When I run univariate models, it works fine and the posteriors
distributions look very good.
Strangely, the more data I have, the more models crash (the largest
amount of data I have is 65000 data for 2400 individuals for one
model).
The model looks like:
priorExp<-list(G=list(G1=list(V=diag(2), nu=2,
alpha.mu=rep(0,2),alpha.V=diag(2)*100000),
G2=list(V=diag(2), nu=2, alpha.mu=rep(0,2),alpha.V=diag(2)*100000),
G3=list(V=diag(2), nu=2, alpha.mu=rep(0,2),alpha.V=diag(2)*100000)),
R=list(V=diag(2), nu=2))
spz<-MCMCglmm(cbind(age1_2,age5_6)~trait-1 + trait:tse+
trait:I(tse^2)+ trait:joe + trait:I(joe^2),
random=~us(trait):animal+us(trait):ID+us(trait):annee ,
rcov=~us(trait):units,nitt=150000, thin=1000, burnin=100000,
prior=priorExp, verbose=TRUE, pedigree=ped,
family=c("gaussian","gaussian"), data=dat)
For the fixed effects : I use 4 continuous parameters as correction
for each trait
For the random effects: I use, individuals, years and animal parameters
I have also tried more informative prior (as described in WAMWIKI)
but the problem was the same.
To give you an example :
Because of computing limitation, I use multi-chain process. I run
several times the same model (as above) and concatenate results
(same prior,same burning, same thin and random seed) to obtain at
least 1000 estimates (50 estimates by model). In this context, I ran
50 bivariable models with the age-class age1_2 and the age-class
age5_6 but only 9 models of the 50 models reached the end.
When we look fixed parameters estimates (estimate are binded for the
nine models : http://ubuntuone.com/3Gi8GwjcRk3P01MxJp2qLe ), we can
see that the estimates are really close to 0. Could it be the problem?
When we look ramdom parameters estimates (estimate are binded for
the nine models : http://ubuntuone.com/42oaP9euG1m2LNipMawcHX ), the
residual estimates do not look very good. Could it be the problem?
Last thing, if I try to add a cubic effect, all the models crash
(same error than before or memory mapped error).
I really do not know where the problem comes from. Do you have an idea?
Thanks
--
Stephane Chantepie
CNRS Phd candidate
Mus?um national d'Histoire naturelle
55 rue Buffon
75005 paris
E-mail : chantepie_at_mnhn.fr
--
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
Dear Jarrod and others, It is not a segfault error, the runs actually stop with one of these two errors. My covariates have already been rescaled (centered to be precise). I have tried to center-reduce the covariates (by usind scale) but the result was the same. Also, when I use cubic term for "tse" and "joe" , the models stop running before 1000 iterations. Could the problem come from the number of fixed parameters to estimate? It seems strange because I have a quiet big data set. kind regards stephane Jarrod Hadfield <j.hadfield at ed.ac.uk> a ?crit?:
Dear Stephane, When you say crash do you mean crash in the sense of a segfault or in the sense that it stops with the errors: -Mixed model equations singular: use a (stronger) prior -ill-conditioned G/R structure: use proper priors if you haven't or rescale data if you have If the latter, it may just require a rescaling of your continuous covariates by using something like scale(). If the former, it would be good for me to have a reproducible example as it means there is a bug. Cheers, Jarrod uoting chantepie at mnhn.fr on Fri, 09 Mar 2012 10:33:42 +0100:
Dear all,
I am using MCMCglmm function to construct bivariate animal models
of bustard sperm production according to age-classes.
The problem is that the models can stochastically crash before the
end of the run (at 2000 iterations or 120000 or other) or can
finish. For the model which does not finish, R returns different
errors as:
-Mixed model equations singular: use a (stronger) prior
-ill-conditioned G/R structure: use proper priors if you haven't or
rescale data if you have
For the models which reach the end, the estimations of genetic
additive variance appear quite good (nice bell shaped posterior
disctribution).
The problem still remains when I remove the animal term.
When I run univariate models, it works fine and the posteriors
distributions look very good.
Strangely, the more data I have, the more models crash (the largest
amount of data I have is 65000 data for 2400 individuals for one
model).
The model looks like:
priorExp<-list(G=list(G1=list(V=diag(2), nu=2,
alpha.mu=rep(0,2),alpha.V=diag(2)*100000),
G2=list(V=diag(2), nu=2, alpha.mu=rep(0,2),alpha.V=diag(2)*100000),
G3=list(V=diag(2), nu=2, alpha.mu=rep(0,2),alpha.V=diag(2)*100000)),
R=list(V=diag(2), nu=2))
spz<-MCMCglmm(cbind(age1_2,age5_6)~trait-1 + trait:tse+
trait:I(tse^2)+ trait:joe + trait:I(joe^2),
random=~us(trait):animal+us(trait):ID+us(trait):annee ,
rcov=~us(trait):units,nitt=150000, thin=1000, burnin=100000,
prior=priorExp, verbose=TRUE, pedigree=ped,
family=c("gaussian","gaussian"), data=dat)
For the fixed effects : I use 4 continuous parameters as correction
for each trait
For the random effects: I use, individuals, years and animal parameters
I have also tried more informative prior (as described in WAMWIKI)
but the problem was the same.
To give you an example :
Because of computing limitation, I use multi-chain process. I run
several times the same model (as above) and concatenate results
(same prior,same burning, same thin and random seed) to obtain at
least 1000 estimates (50 estimates by model). In this context, I
ran 50 bivariable models with the age-class age1_2 and the
age-class age5_6 but only 9 models of the 50 models reached the end.
When we look fixed parameters estimates (estimate are binded for
the nine models : http://ubuntuone.com/3Gi8GwjcRk3P01MxJp2qLe ), we
can see that the estimates are really close to 0. Could it be the
problem?
When we look ramdom parameters estimates (estimate are binded for
the nine models : http://ubuntuone.com/42oaP9euG1m2LNipMawcHX ),
the residual estimates do not look very good. Could it be the
problem?
Last thing, if I try to add a cubic effect, all the models crash
(same error than before or memory mapped error).
I really do not know where the problem comes from. Do you have an idea?
Thanks
--
Stephane Chantepie
CNRS Phd candidate
Mus?um national d'Histoire naturelle
55 rue Buffon
75005 paris
E-mail : chantepie_at_mnhn.fr
--
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
1 day later
Hi Stephane, Its very hard to know what could be causing the problem with out having the data. The three things that you could try are: a) scale the response variable b) work out if some terms really are weakly identified (MCMCglmm drops non-identified terms in the fixed effects by default) c) the error can appear I think if the residual covariance matrix becomes singular. This could happen, amongst other reasons, because non-modelled non-genetic effects are contributing to the resemblance between relatives and inflating Va by so much that Ve goes to zero. Either i) running it for some iterations before it terminates or ii) fit the model in ASReml should give you some idea which of these problems are most likely. Cheers, Jarrod Quoting chantepie at mnhn.fr on Fri, 09 Mar 2012 17:33:52 +0100:
Dear Jarrod and others, It is not a segfault error, the runs actually stop with one of these two errors. My covariates have already been rescaled (centered to be precise). I have tried to center-reduce the covariates (by usind scale) but the result was the same. Also, when I use cubic term for "tse" and "joe" , the models stop running before 1000 iterations. Could the problem come from the number of fixed parameters to estimate? It seems strange because I have a quiet big data set. kind regards stephane Jarrod Hadfield <j.hadfield at ed.ac.uk> a ?crit?:
Dear Stephane, When you say crash do you mean crash in the sense of a segfault or in the sense that it stops with the errors: -Mixed model equations singular: use a (stronger) prior -ill-conditioned G/R structure: use proper priors if you haven't or rescale data if you have If the latter, it may just require a rescaling of your continuous covariates by using something like scale(). If the former, it would be good for me to have a reproducible example as it means there is a bug. Cheers, Jarrod uoting chantepie at mnhn.fr on Fri, 09 Mar 2012 10:33:42 +0100:
Dear all,
I am using MCMCglmm function to construct bivariate animal models
of bustard sperm production according to age-classes.
The problem is that the models can stochastically crash before the
end of the run (at 2000 iterations or 120000 or other) or can
finish. For the model which does not finish, R returns different
errors as:
-Mixed model equations singular: use a (stronger) prior
-ill-conditioned G/R structure: use proper priors if you haven't
or rescale data if you have
For the models which reach the end, the estimations of genetic
additive variance appear quite good (nice bell shaped posterior
disctribution).
The problem still remains when I remove the animal term.
When I run univariate models, it works fine and the posteriors
distributions look very good.
Strangely, the more data I have, the more models crash (the
largest amount of data I have is 65000 data for 2400 individuals
for one model).
The model looks like:
priorExp<-list(G=list(G1=list(V=diag(2), nu=2,
alpha.mu=rep(0,2),alpha.V=diag(2)*100000),
G2=list(V=diag(2), nu=2, alpha.mu=rep(0,2),alpha.V=diag(2)*100000),
G3=list(V=diag(2), nu=2, alpha.mu=rep(0,2),alpha.V=diag(2)*100000)),
R=list(V=diag(2), nu=2))
spz<-MCMCglmm(cbind(age1_2,age5_6)~trait-1 + trait:tse+
trait:I(tse^2)+ trait:joe + trait:I(joe^2),
random=~us(trait):animal+us(trait):ID+us(trait):annee ,
rcov=~us(trait):units,nitt=150000, thin=1000, burnin=100000,
prior=priorExp, verbose=TRUE, pedigree=ped,
family=c("gaussian","gaussian"), data=dat)
For the fixed effects : I use 4 continuous parameters as
correction for each trait
For the random effects: I use, individuals, years and animal parameters
I have also tried more informative prior (as described in WAMWIKI)
but the problem was the same.
To give you an example :
Because of computing limitation, I use multi-chain process. I run
several times the same model (as above) and concatenate results
(same prior,same burning, same thin and random seed) to obtain at
least 1000 estimates (50 estimates by model). In this context, I
ran 50 bivariable models with the age-class age1_2 and the
age-class age5_6 but only 9 models of the 50 models reached the end.
When we look fixed parameters estimates (estimate are binded for
the nine models : http://ubuntuone.com/3Gi8GwjcRk3P01MxJp2qLe ),
we can see that the estimates are really close to 0. Could it be
the problem?
When we look ramdom parameters estimates (estimate are binded for
the nine models : http://ubuntuone.com/42oaP9euG1m2LNipMawcHX ),
the residual estimates do not look very good. Could it be the
problem?
Last thing, if I try to add a cubic effect, all the models crash
(same error than before or memory mapped error).
I really do not know where the problem comes from. Do you have an idea?
Thanks
--
Stephane Chantepie
CNRS Phd candidate
Mus?um national d'Histoire naturelle
55 rue Buffon
75005 paris
E-mail : chantepie_at_mnhn.fr
--
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
4 days later
Dear all, I have tried the different solutions you proposed. I have tried to scale fixed parameters and the response variables between 0 and 1 but the problem remained However, using ASREML, it seems that the problem comes from the residual covariance matrix. As supposed by Jarrod the residual covariance matrix is singular. I am pretty sure that it is due the structure of my data but I really don't know how I can fix this problem. I use several age classes (which represent my traits) and I have intra-annual repeated measurements for several years and for each individual. The fixed parameters I use are: -tse : time since the last ejaculation collect: it is used to take into account the pressure due to the repeated collect. -joe : day of ejaculation : represents the time since the first ejaculation of the year. I use this parameter to take into account the seasonal variation of sperm production. The model is AgeClass1 AgeClass2 ~ at.level(trait,1): tse+ at.level(trait,2): tse+ at.level(trait,1): joe + at.level(trait,2): joe, random=~us(trait):animal+us(trait):ID+us(trait):Year, rcov = ~us(trait):units With the data structure I have, it is impossible to have measurements of the same individual on the same line (snapshot to help comprehension: http://ubuntuone.com/1W19vErC5jUtHqMn1dtmPg ), likely making it impossible to estimate a residual covariance. One of the issues is that I can not see really how to change the structure of the data so it?s estimable. Is there a solution? I?m afraid if I fix the residuals covariance matrix to 0, I will inflate Va. Do you have an idea? many thanks for your help stephane