Re: [R-sig-phylo] How to incorporate intraspecific variation in MCMCglmm
cc-ed back to the list
Some observations seem to have no sampling variation and so have a
residual variance of zero - which is not allowed, and doesn't make much
sense.
Cheers,
Jarrod
On 16/07/2017 00:16, Diogo B. Provete wrote:
Hi Jarrod,
I have tried all the suggestions for specifying the random effects.
However, your last suggestion of making the residual covariance seems
more appropriated in my case, since I have the standard error
calculated already (in terms of jumped distance) and it's interpreted
as intraspecific variance due to sampling error. But when I make my
model as:
model1<-MCMCglmm(mean_distance~type_arena*type_of_stimulus,
*random=~Species, rcov=~idh(units):units*, data=df_spe,
family="gaussian", ginverse = list(Especie=treeAinv), nodes="ALL",
prior=ve_priors, nitt=30, burnin=25000, thin = 100, verbose=FALSE)
with priors:
ve_priors = list(R = list(V = diag(df_spe$se_distance^2), fix=1),
G=list(G1=list(V=1, nu=0.02)))
I get the following error message:
Error in priorformat(if (NOpriorG) { :
V is not positive definite for some prior$G/prior$R elements .
I tried to google it and looked at your course notes ch. 3 and 4, but
couldn't find anything helpful to understand it.
I'm attaching the first few lines of the prior V matrix.
Thanks once again,
Diogo
Em sex, 14 de jul de 2017 às 18:47, Diogo B. Provete
> escreveu:
Dear Jarrod,
Thanks very much for the quick reply.
I'll try to implement the changes in the model.
Have a nice weekend,
Diogo
Em Sex, 14 de jul de 2017 17:48, Jarrod Hadfield
> escreveu:
Hi Diogo,
First, your model1 is unlikely to be valid unless the residual
variance
happens to be 1. You should not fix it at one, and use a prior
like:
prior = list(R = list(V = 1, nu = 0.02), G=list(G1=list(V=1,
nu=0.02)))
Note that the residual variance (Ve) is the intra-specific
variance,
assumed constant over species, as in Lynch's h2/Pagels Lambda.
If you
have the standard error (se) of each observation you can also
allow
heterogeneity between observations as in random-effect
meta-analysis:
random=~Specie+idh(se)
If you fix the variance associated with se at 1 in the prior
prior = list(R = list(V = 1, nu = 0.02), G=list(G1=list(V=1,
nu=0.02), G2=list(V=1, fix=1)))
Then the intraspecific variance for species i is se^2_i+Ve
If you do not fix the variance associated with se at 1, for
example with
prior = list(R = list(V = 1, nu = 0.02), G=list(G1=list(V=1,
nu=0.02), G2=list(V=1, nu=0.02)))
and estimate the associated variance (Vse). Then the
intraspecific variance
for species i is Vse*se^2_i+Ve. This is useful if you only
know the standard error up to proportionality.
Alternatively you can fit fixed-effect meta-analysis where all
the intraspecific variance is presumed to be due to sampling
error:
random=~Specie, rcov=~idh(units):units
with prior:
prior = list(R = list(V = diag(se^2), fix=1),
G=list(G1=list(V=1, nu=0.02)))
The intraspecific variance for species i is se^2_i. This is
quite an inefficient way of doing fixed-effect meta-analysis,
and one day I will make such analyses easier to fit, and
quicker to fit..
Also, don't use nodes="TIPS" in the call to inverseA. I only
allowed this option because its the parameterisation used by
most other software, but its a really bad way of doing it.
Stick with the default, node="ALL".
Cheers,
Jarrod
On 14/07/2017 14:26, Diogo B. Provete wrote:
> treeAinv<-inverseA(phylo,nodes="TIPS",scale=TRUE)$Ainv
>
> I included the following priors:
>
> prior = list(R = list(V = 1, fix = 1), G=list(G1=list(V=1,
nu=0.02)))
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
--
*Diogo B. Provete, PhD.*
FAPESP Post-doctoral fellow
Department of Environmental Sciences
Centre for Sciences and Technologies for Sustainability
Federal University of São Carlos
Sorocaba
Gothenburg Global Biodiversity Centre
Box 462 SE-405 30
Göteborg, Sverige
diogoprovete.weebly.com
Cell phone: +5515981022137
Skype: diogoprovete
Editor: Amphibia-Reptilia | Biodiversity Data Journal
--
*Diogo B. Provete, PhD.*
FAPESP Post-doctoral fellow
Department of Environmental Sciences
Centre for Sciences and Technologies for Sustainability
Federal University of São
Re: [R-sig-phylo] How to incorporate intraspecific variation in MCMCglmm
Hi Diogo, First, your model1 is unlikely to be valid unless the residual variance happens to be 1. You should not fix it at one, and use a prior like: prior = list(R = list(V = 1, nu = 0.02), G=list(G1=list(V=1, nu=0.02))) Note that the residual variance (Ve) is the intra-specific variance, assumed constant over species, as in Lynch's h2/Pagels Lambda. If you have the standard error (se) of each observation you can also allow heterogeneity between observations as in random-effect meta-analysis: random=~Specie+idh(se) If you fix the variance associated with se at 1 in the prior prior = list(R = list(V = 1, nu = 0.02), G=list(G1=list(V=1, nu=0.02), G2=list(V=1, fix=1))) Then the intraspecific variance for species i is se^2_i+Ve If you do not fix the variance associated with se at 1, for example with prior = list(R = list(V = 1, nu = 0.02), G=list(G1=list(V=1, nu=0.02), G2=list(V=1, nu=0.02))) and estimate the associated variance (Vse). Then the intraspecific variance for species i is Vse*se^2_i+Ve. This is useful if you only know the standard error up to proportionality. Alternatively you can fit fixed-effect meta-analysis where all the intraspecific variance is presumed to be due to sampling error: random=~Specie, rcov=~idh(units):units with prior: prior = list(R = list(V = diag(se^2), fix=1), G=list(G1=list(V=1, nu=0.02))) The intraspecific variance for species i is se^2_i. This is quite an inefficient way of doing fixed-effect meta-analysis, and one day I will make such analyses easier to fit, and quicker to fit.. Also, don't use nodes="TIPS" in the call to inverseA. I only allowed this option because its the parameterisation used by most other software, but its a really bad way of doing it. Stick with the default, node="ALL". Cheers, Jarrod On 14/07/2017 14:26, Diogo B. Provete wrote: treeAinv<-inverseA(phylo,nodes="TIPS",scale=TRUE)$Ainv I included the following priors: prior = list(R = list(V = 1, fix = 1), G=list(G1=list(V=1, nu=0.02))) -- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] How to incorporate intraspecific variation in MCMCglmm
Dear list, I'm trying to fit a mixed-effects model to a problem that includes: a continuous response variable (mean jumped distance by 17 species of frogs) in 3 different arenas and under 2 types of stimuli. So, it's a kind of two-way ANOVA design. I have the species phylogeny and prepared its inverse: treeAinv<-inverseA(phylo,nodes="TIPS",scale=TRUE)$Ainv I included the following priors: prior = list(R = list(V = 1, fix = 1), G=list(G1=list(V=1, nu=0.02))) and the model: model1<-MCMCglmm(mean_distance~type_arena*type_of_stimulus, random=~Specie, data=df_spe, family="gaussian", ginverse = list(Specie=treeAinv), nodes="ALL", prior=prior, nitt=30, burnin=25000, thin = 100, verbose=FALSE) which seems to converge well (heidel.diag p-value = 0.428), ESS greater than 1000 for both random and fixed effects. My question is: how can I incorporate intraspecific variation in the jumped distance? I have the standard error of the distance and tried to make the model as: random=~Species+ us(1+se_distance):Species but I keep getting an error about prior specification: Error in MCMCglmm(mean_distance ~ type_arena * type_og_stimulus, random = ~Specie + : prior$G has the wrong number of structures Thank you in advance, Diogo -- *Diogo B. Provete, PhD.* FAPESP Post-doctoral fellow Department of Environmental Sciences Centre for Sciences and Technologies for Sustainability Federal University of São Carlos Sorocaba Gothenburg Global Biodiversity Centre Box 462 SE-405 30 Göteborg, Sverige diogoprovete.weebly.com Cell phone: +5515981022137 Skype: diogoprovete Editor: Amphibia-Reptilia | Biodiversity Data Journal [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
