Hi,
The least parsimonious model (and not the one I would necessarily
recommend fitting) is:
m1<-MCMCglmm(cbind(X1,X2,X3,Y)~trait,
random=~us(trait):species+us(trait):species.ide,
rcov=~us(trait):units,
ginverse=list(species=tree))
where species and species.ide are columns of species names.
This deals with the measurement error on the species means, and also
allows you to address the fact that the regressions of the X's on Y may be
different at different levels. The method advocated by van de Pol has the
problem that the mean in the mean centering is just the observed mean
rather than the true unobserved mean. For example, imagine that you only
had one observation for some of the species. You can obtain the regression
coefficients at each level, by using the relationship beta =
VAR(X)^{-1}COV(X,Y). For example, the posterior distribution of the
regression coefficients at the phylogenetic level would be:
reg.coef<-function(x, X=1:3, Y=4){
V<-matrix(x,c(X,Y),c(X,Y))
solve(V[X,X], V[X,Y])
}
apply(m1$VCV[,1:16], 1, reg.coef)
The model doesn't deal with measurement error on the individual
measurements, but if you had repeat measurements per individual you could
also fit these (as a diagonal matrix, rather than unstructured).
After taking into account measurement error, some people suggest that
species.ide should be dropped from the model. I am not completely convinced
by this argument.
Priors are going to be a pain in this model.
You could replace the us structures by ante3 structures. The model is then
fitted directly in terms of the regression coefficients. Antedependence
regression coefficients 3,5,6 are the regressions of X3, X2 and X1 on Y. If
you are interested in this we have a mini-tutorial associated with a
recently submitted paper I can send you.
Cheers,
Jarrod
Quoting Alberto Gallano <alberto.gc8 at gmail.com> on Sun, 3 Jan 2016
15:45:26 -0500:
Hi Jarrod,
yes, that's right, I have multiple measurements for both response and
predictors and these are measured on the same individuals. The model i'm
fitting is very similar to the model called "model_repeat2" from Modern
Phylogenetic Comparative Methods:
http://www.mpcm-evolution.org/practice/online-practical-material-chapter-11/chapter-11-2-multiple-measurements-model-mcmcglmm
same random effects structure, same between/within structure for the fixed
effects, and i'm also using the inverse of the matrix of phylogenetic
correlation.
@Dimitri: I'm aware of the de Villemereuil et al. approach, which, If I
understand correctly, does a version of orthogonal regression (in JAGS).
I'm trying find out if this is possible in MCMCglmm.
best,
Alberto
On Sun, Jan 3, 2016 at 12:05 PM, Dimitri Skandalis <
da.skandalis at gmail.com>
wrote:
Hi Alberto,
When you say you have multiple observations for each species, do you
mean
that you have multiple observations for the response and the
predictors? Do
you expect the response and/or the predictors to be correlated at the
observation level (for example are they measured on the same
individuals)?
I presume the answer to both these questions is yes if you wish to use
the
van de Pol method?
Cheers,
Jarrod
Quoting Alberto Gallano <alberto.gc8 at gmail.com> on Sun, 3 Jan 2016
10:35:02 -0500:
Hi Jarrod,
I don't know the measurement error in the predictors in advance, so I
guess
it would need to be estimated simultaneously. I'm not 100% sure what
you
mean by 'multiple observations for each predictor variable'. I have
data
on
132 species and have multiple observations (7 to 80) for each species.
I'm
using a species level random effect and a phylogenetic covariance
matrix
(using ginverse) to account for phylogenetic autocorrelation, and I'm
also
using van de Pol and Wright's (2009) method for partitioning slopes
into
between- and within-species (i'm interested in the between species
slope).
My understanding is that neither of these things fits a model in which
orthogonal residuals are minimized.
best,
Alberto
On Sun, Jan 3, 2016 at 5:24 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk>
wrote:
Hi Alberto,
Do you know the measurement error in the predictors in advance or do
you
have multiple observations for each predictor variable and wish to
estimate
the error simultaneously?
Cheers,
Jarrod
Quoting Malcolm Fairbrother <M.Fairbrother at bristol.ac.uk> on Sat, 2
Jan
2016 14:47:08 -0800:
Dear Alberto (I believe),
To my knowledge, this is not possible in MCMCglmm (though Jarrod
Hadfield,
the package author, may weigh in with another response).
A collaborator and I have been working on a paper that shows how to
fit
such models in JAGS (and perhaps Stan), though thus far we've only
been
able to fit such models correcting for measurement error in the
predictors
at the lowest level. Multiple such predictors (including with
different
measurement error variances) are no problem.
That paper, however, is probably still some months away from being
finished
and presentable. In the meantime, I don't know of any good options
for
you.
If other subscribers to this list have any ideas, I'll be quite
interested
too!
- Malcolm
Date: Tue, 29 Dec 2015 16:09:53 -0500
From: Alberto Gallano <alberto.gc8 at gmail.com>
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] MCMCglmm error-in-variables (total least
squares)
model?
I posted this question on Stack Overflow a week ago but received no
answers:
http://stackoverflow.com/questions/34446618/bayesian-error-in-variables-total-least-squares-model-in-r-using-mcmcglmm
This may be a more appropriate venue.
I am fitting some Bayesian linear mixed models using the MCMCglmm
package.
My data includes predictors that are measured with error. I'd
therefore
like to build a model that takes this into account. My understanding
is
that a basic mixed effects model in MCMCglmm will minimize error
only
for
the response variable (as in frequentist OLS regression). In other
words,
vertical errors will be minimized. Instead, I'd like to minimize
errors
orthogonal to the regression line/plane/hyperplane.
1. Is it possible to fit an error-in-variables (aka total least
squares)
model using MCMCglmm or would I have to use JAGS / STAN to do
this?
2. Is it possible to do this with multiple predictors in the same
model
(I have some models with 3 or 4 predictors, each measured with
error)?
[[alternative HTML version deleted]]