[R-sig-ME] lme with a known correlation matrix for the random effects

[Vmitov@gmail.com]

Dear Thierry,

Thanks a lot for letting me know this!

Venelin

> On 12.12.2014, at 10:54, ONKELINX, Thierry <Thierry.ONKELINX at inbo.be> wrote:
> 
> Dear Mitov,
> 
> The nlme package allows only for correlated residuals within the lowest level of random effects. If you have the patient id as random effect, then only the correlation of the residuals within each patient can be modeled (or set to a predefined value with fixed = TRUE). Residuals among patients are assumed to be independent.
> 
> The INLA packages allows for correlated random effects.
> 
> Best regards
> 
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
> Kliniekstraat 25
> 1070 Anderlecht
> Belgium
> + 32 2 525 02 51
> + 32 54 43 61 85
> Thierry.Onkelinx at inbo.be
> www.inbo.be
> 
> To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
> ~ Sir Ronald Aylmer Fisher
> 
> The plural of anecdote is not data.
> ~ Roger Brinner
> 
> The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
> ~ John Tukey
> 
> 
> -----Oorspronkelijk bericht-----
> Van: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-project.org] Namens Mitov Venelin
> Verzonden: donderdag 11 december 2014 22:26
> Aan: r-sig-mixed-models at r-project.org
> Onderwerp: [R-sig-ME] lme with a known correlation matrix for the random effects
> 
> Dear R-sig mailing list,
> 
> Currently, I'm trying to use the function lme from the nlme package to model epidemiological data. My question seems to be a specific one and I couldn't find any similar model formulation using the lme function. That's why I'm writing you.
> 
> The data I have is as follows:
> 
> 1. y: a numeric response vector consisting of one phenotype measurement for each of N patients.
> 
> 2. A known known correlation matrix, Psi, between the patients.
> 
> I would like to model this data as a linear mixed effects model as follows:
> 
> y_i = beta + g_i + e_i,
> 
> i=1,...,N
> 
> Here are the model assumptions:
> 
> beta is a fixed population grand-mean (i.e. intercept);
> 
> g_i is the pathogen-specific contribution to y_i and is considered as a random effect that is correlated along patients (explained in more detail below);
> 
> e_i is a patient-specific random effect including the effects of patient's immune system, any possible interaction between the patient's immune system and the pathogen, and measurement error for y_i.
> 
> In matrix-notation the model can be written as follows
> 
> y = 1*beta + I_NxN*g + I_NxN*e
> 
> In the formula above 1 denotes a N-dimensional vector of ones, I_NxN denotes the NxN identity matrix.
> The random N-dimensional vector g is assumed to be normally distributed with mean 0 and a known correlation matrix, Psi.
> The random N-dimensional vector e is assumed to be normally distributed with mean 0 and and identity correlation matrix I_NxN.
> 
> The model has three parameters:
> - beta;
> - sigma_g, such that covariance matrix of g is equal to Psi*sigma_g^2 (Psi is known and fixed);
> - sigma_e such that the covariance matrix of e is I_NxN*sigma_e^2;
> 
> Note that in the lme terminology each of the patients represents a group by itself and there is 1 observation per group.
> I'm able to fit the model using maximum likelihood, but I'm trying to use the lme function mostly because of its REML fitting method.
> 
> So far, I've tried to use the lme function but I'm getting errors and I don't know how to formulate my model in the lme syntax.
> 
> N=5
> y <- c(7.38, 7.25, 7.34, 7.30, 7.06)
> names(y) <- c("t4", "t2", "t5", "t1", "t3") g <- factor(names(y)) data <- data.frame(y=y, g=g)
> 
> Psi <- matrix(c(1, 0.00, 0.00, 0.00, 0.00,
>                0, 1.00, 0.28, 0.25, 0.78,
>                0, 0.28, 1.00, 0.83, 0.26,
>                0, 0.25, 0.83, 1.00, 0.23,
>                0, 0.78, 0.26, 0.23, 1.00), nrow=N, ncol=N)
> rownames(Psi)<-colnames(Psi)<-names(y)
> 
> pdMatPsi <- pdMat(Psi, form=~1|g, pdClass='pdSymm') # Not sure if the formula is correct and if pdSymm is the right pdMat class to use.
> 
> 
>> lme(y~1, data, random=pdMatPsi)
> Error in getGroups.data.frame(dataMix, groups) :
>  invalid formula for groups
> 
> For this error I guess the problem is that lme tries to look for a formula in the data.frame data which, as I've read in a vignette, is normally used only for plotting the data, so I'd rather keep the data.frame and not extend it to a groupedData.
> 
>> lme(y~1, data, random=list(g=pdMatPsi))
> Error in matrix(unlist(value), nrow = nrow(data), dimnames = list(row.names(data),  :
>  length of 'dimnames' [2] not equal to array extent In addition: Warning message:
> In Ops.factor(1, g) : | not meaningful for factors
> 
> Here I don't understand what am I doing wrong. I've tried other formulas in the call to pdMat but I'm still getting errors.
> Is it possible to use lme to fit the above model to the data? Also, is it possible to use the more recent lme4 package as it is supposed to be much faster?
> 
> Thank you and best regards,
> MV
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> Disclaimer Bezoek onze website / Visit our website<https://drupal.inbo.be/nl/disclaimer-mailberichten-van-het-inbo>