Skip to content

Adding additional factor to previously working MCMCglmm model kills it. Problem with priors? Help?

3 messages · Jarrod Hadfield, Tricia Markle

Original
Tricia Markle
# 
Hello,



I have a working MCMCglmm model with phylogenetic consideration and repeat
measures. I realized after the fact that ?species? wasn?t properly included
in the model. When I added this additional factor (of 17 levels), however,
I received the following error message:

--------



Error in MCMCglmm(LVO2 ~ 1 + Temp + Acclm + Lat_Ext + LMass + Sex + species
+  :  ill-conditioned G/R structure: use proper priors if you haven't or
rescale data if you have

In addition: Warning message:

In MCMCglmm(LVO2 ~ 1 + Temp + Acclm + Lat_Ext + LMass + Sex + species +  :

  some fixed effects are not estimable and have been removed. Use
singular.ok=TRUE to sample these effects, but use an informative prior!



----------



I had some help with my priors originally so I am not sure the best way to
now tweak them to make them work? Could different priors help or is there
something ?wrong? with adding species as another factor?



My hypothesis centers around the remaining significant interaction term
Acclm:Lat_Ext ? whether the relationship of acclimation on oxygen
consumption (VO2) is differs depending on the latitudinal extent of a
species (while also considering a number of covariates).



Here is my code:



library(ape)

library(MCMCglmm)

dataset<-read.csv(file="RespData.csv", head=TRUE)

attach(dataset)

str(dataset) # confirming that sex, range, species, and ID are all factors



#Phylogeny Component

tree<-read.tree("Plethodontidae_comb61_PL.phy")

species<-c("D._carolinensis_KHK103", "D._fuscus_KHK142",
"D._ochrophaeus_WKS05", "D._ocoee_B_KHK62", "D._orestes_KHK129",
"D._monticola_A",  "D._santeetlah_11775", "P_cinereus", "P_cylindraceus",
"P_glutinosus", "P_hubrichti", "P_montanus", "P_punctatus", "P_richmondi",
"P_teyahalee", "P_virginia", "P_wehrlei")

pruned.tree<-drop.tip(tree,tree$tip.label[-match(species,
tree$tip.label)])# Prune tree to just include species of interest

sptree<-makeNodeLabel(pruned.tree, method="number", prefix="node")



treeAinv<-inverseA(sptree, nodes="TIPS")$Ainv

random=~us(1+Temp):species

prior<-list(G=list(G1=list(V=diag(2), nu=2, alpha.mu=c(0,0),
alpha.V=diag(2)*1000)), R=list(V=diag(1), nu=0.002))



#Model 1

model1<MCMCglmm(LVO2~1+Temp+Acclm+Lat_Ext+LMass+Sex+species+Temp*Acclm+Temp*Lat_Ext+Acclm*Lat_Ext+Temp*Acclm*Lat*Ext,
random=random, data=dataset, family="gaussian",
ginverse=list(species=treeAinv), prior=prior, nitt=300000, burnin=25000,
thin = 1000, verbose=FALSE)



Thank you kindly for your help!



Tricia
Jarrod Hadfield
# 
Hi,

If you want iid species effects then put them in the random part. You  
will need to rename them so that they are not associated with the  
ginverse. For example,

dataset$species.ide<-dataset$species

and fit species.ide as a random effect.

You only have 17 species. This is not enough to get precise estimates  
of the variance components and so you should expect prior sensitivity,  
particularly with respect to the phylogenetic part. I would try and  
simplify your model. I'm not sure how many of the fixed effects are  
species-level or individual-level but I would try and reduce the  
complexity of the fixed part of the model too (for example is the  
4-way interaction needed?). The warning usually indicates that the  
model is overparameterised.

Cheers,

Jarrod




  Quoting Tricia Markle <markl033 at umn.edu> on Wed, 29 Apr 2015 00:24:53 -0500:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Tricia Markle
    

    

      
      #
    

  


  

      
Hi Jarrod,

Thanks for the response.

If I am understanding correctly, would I then end up with the code below?
If species is incorporated in the random part, would I then not include it
as a fixed factor?

tree<-read.tree("Plethodontidae_comb61_PL.phy")
species<-c("D._carolinensis_KHK103", "D._fuscus_KHK142",
"D._ochrophaeus_WKS05", "D._ocoee_B_KHK62", "D._orestes_KHK129",
"D._monticola_A",  "D._santeetlah_11775", "P_cinereus", "P_cylindraceus",
"P_glutinosus", "P_hubrichti", "P_montanus", "P_punctatus", "P_richmondi",
"P_teyahalee", "P_virginia", "P_wehrlei")
pruned.tree<-drop.tip(tree,tree$tip.label[-match(species,
tree$tip.label)])# Prune tree to just include species of interest
sptree<-makeNodeLabel(pruned.tree, method="number", prefix="node") #rename
nodes to be unique
treeAinv<-inverseA(sptree, nodes="TIPS")$Ainv

dataset$species.ide<-dataset$species
random=~us(1+Temp):species.ide

prior<-list(G=list(G1=list(V=diag(2), nu=2, alpha.mu=c(0,0),
alpha.V=diag(2)*1000)), R=list(V=diag(1), nu=0.002))

model6<-MCMCglmm(LVO2~1+Temp+Acclm+Range+Mass+Sex+Acclm*Range,
random=random, data=dataset, family="gaussian",
ginverse=list(species=treeAinv), prior=prior, nitt=300000, burnin=25000,
thin = 1000, verbose=FALSE)

Alternatively, someone mentioned to me that to reduce model complexity I
could group the handful of species I have with smaller samples sizes and
label them all as "other" in a second species list ("New_spp"). When I
added this new species assortment as a fixed factor in my original model it
runs! Is this a reasonable thing to do?

Thanks, Tricia

On Wed, Apr 29, 2015 at 12:46 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk>
wrote: