[R-sig-ME] weighting nlme and multivariate outcomes

Thu Aug 20 22:12:35 CEST 2009

Dear Fabian and interested others,

I thought some more about the use of corSymm to accomplish your goal, as this is a question I was also thinking about myself awhile back.  The model below is close but one problem still needs to be worked out it seems.  Any other input much appreciated.  Summary below:

Consider the case where we have n1 subjects and two variables U and V, each measured at n2 time points.   Assume each variable arises from a random-intercept, random-slope model, with all random effects correlated.  Thus, the random-effects covariance matrix is 4x4.

Regarding the residual error term, for illustration suppose there are n2=2 time points.  Then for a given subject the response vector Y_i = (U_i', V_i')'  and residual error vector (e_i, d_i')' are both 4x1.  The covariance matrix of the residual error vector will thus be 4x4: the upper 2x2 being the covariance matrix of e_i, the lower 2x2 being the covariance matrix of d_i, and the off-diagonal 2x2 being the covariance matrix of e_i and d_i.  For instance, the two diagonal elements of this last 2x2 represent the covariance of the error terms of the different variables at time 1 and time 2 (the corresponding off-diagonal elements represent the covariance when each variable is measured at a _different_ time point).

Now, here are various ways to fit models to these data under different assumptions for the random-effects; the residual errors are fit with different variances per variable, and are initially assumed independent within and between variables (first copy and past the code block at the end of this e-mail that simulates the data) .

# all random-effects correlated; different residuial error variance per variable
mv1.lme = lme(y ~ -1 + U.ind + U.time + V.ind + V.time, random = ~ -1 + U.ind + U.time + V.ind + V.time | subject,
            data = data.mv.grp, control=nlmeControl(msMaxIter = 1000), weights = varIdent(form = ~ 1 | U.ind))
# all random-effects independent; need to have data as grouped-data object for pdDiag to work below
mv2.lme = lme(y ~ -1 + U.ind + U.time + V.ind + V.time, random = pdDiag(~ -1 +  U.ind + V.ind + V.time + U.time  ),
            data = data.mv.grp, control=nlmeControl(msMaxIter = 1000), weights = varIdent(form = ~ 1 | U.ind))
# random effects not correlated across variables
mv3.lme = update(mv2.lme, random = list(subject = pdBlocked(list(pdSymm(~ -1 +  U.ind + U.time),
             pdSymm(~ -1 + V.ind + V.time)))))

To relax the assumption of independence of residual errors across variables, it seems that corSymm should be sufficient.  Below I initialize the values for corSymm to correspond to the values used in the simulation, which are:

> Sig.corr
     [,1] [,2] [,3] [,4]
[1,]  1.0  0.0  0.1  0.0
[2,]  0.0  1.0  0.0  0.1
[3,]  0.1  0.0  1.0  0.0
[4,]  0.0  0.1  0.0  1.0

Thus,

value.star = offdiag.to.vector(Sig.corr)[[1]]  ## use initial values as actual values form simulation
cs1Symm = corSymm(value = value.star, form = ~ 1 | subject)
cs1Symm = Initialize(cs1Symm, data = data.mv)

However, regardless of which model above is used, adding correlation = cs1Symm introduces convergence problems:

mv3.lme = update(mv1.lme, correlation = cs1Symm, control=nlmeControl(msMaxIter = 1000))
Error in lme.formula(fixed = y ~ -1 + U.ind + U.time + V.ind + V.time,  :
  nlminb problem, convergence error code = 1
  message = false convergence (8)

It would seem that we need to perhaps fix some of the values of the correlation matrix to 0, viz., all off-diagonal elements of each of the 2x2 blocks mentioned earlier, but I do not see a way to do that based on the description in P&B (p.235).  It seems like this exists for varFunc but not for corStruct.  Does anyone know a way around this to reduce the number of parameters that are being estimated for the residual error accordingly?

Perhaps the design values of my simulation are part of the problem but I played around with them and didn't have much success.  If this last issue with corSymm can be solved it seems that a wide class of multivariate models will be able to be fit with lme.

Cheers,
David

##############

library("MASS"); library("nlme"); library("lattice")
set.seed(100); n1 = 50; n2 = 4
# U: random-intercept b1, random-slope b2
sigma.b1.sq = 3; sigma.b2.sq = 8;
# V: random-intercept c1, random-slope c2
sigma.c1.sq = 6; sigma.c2.sq = 12;
# correlations between random-effects
corr.b1.c1 = .1; sigma.b1.c1 = corr.b1.c1*(sqrt(sigma.b1.sq)*(sqrt(sigma.c1.sq)))
corr.b1.b2 = .2; sigma.b1.b2 = corr.b1.b2*(sqrt(sigma.b1.sq)*(sqrt(sigma.b2.sq)))
corr.c1.c2 = .2; sigma.c1.c2 = corr.c1.c2*(sqrt(sigma.c1.sq)*(sqrt(sigma.c2.sq)))
corr.b1.c2 = .05; sigma.b1.c2 = corr.b1.c2*(sqrt(sigma.b1.sq)*(sqrt(sigma.c2.sq)))
corr.b2.c1 = .05; sigma.b2.c1 = corr.b2.c1*(sqrt(sigma.b2.sq)*(sqrt(sigma.c1.sq)))
corr.b2.c2 = .35; sigma.b2.c2 = corr.b2.c1*(sqrt(sigma.b2.sq)*(sqrt(sigma.c2.sq)))
# residual error term:
sigma.u.sq = 14; sigma.v.sq = 25
corr.u.v = .1; sigma.u.v = corr.u.v*sqrt(sigma.u.sq)*sqrt(sigma.v.sq)
# simulate multivariate longitudinal profile
U = V = matrix(0, n2, n1)
for (i in 1:n1) {
    b.i = mvrnorm(n = 1, mu = c(0,0,0,0), Sigma = matrix(c(
            sigma.b1.sq, sigma.b1.b2, sigma.b1.c1,  sigma.b1.c2,
            sigma.b1.b2, sigma.b2.sq, sigma.b2.c1,  sigma.b2.c2,
            sigma.b1.c1, sigma.b2.c1, sigma.c1.sq,  sigma.c1.c2,
            sigma.b1.c2, sigma.b2.c2, sigma.c1.c2,  sigma.c2.sq), nrow = 4, ncol = 4)) # correlated random effects
    for (j in 1:n2) {
        eps =  mvrnorm(n = 1, mu = c(0,0), Sigma = matrix(c(sigma.u.sq, sigma.u.v, sigma.u.v, sigma.v.sq), nrow=2,ncol=2))
        U[j,i] = 75 + (13 + b.i[2]) * j + b.i[1] + eps[1]
        V[j,i] = 50 + (2 + b.i[4]) * j + b.i[3] + eps[2]  ## fixed intercept and sloped effects + random-interecept + random-slope
        }
    }
data.mv = data.frame(subject = rep(seq(1,n1), each = n2, 2), y = c(as.vector(U), as.vector(V)),
            U.ind = rep( c(1,0), each = n1*n2), V.ind = rep( c(0,1), each = n1*n2), U.time = c(rep(seq(1:n2), n1), rep(0, n1*n2)), V.time = c(rep(0, n1*n2), rep(seq(1,n2), n1)))
data.mv$time = data.mv$U.time + data.mv$V.time
data.mv.grp = groupedData(y ~ time | subject, data = data.mv)

Sig.upper.left = sigma.u.sq*diag(n2); Sig.lower.right = sigma.v.sq*diag(n2)
Sig.off.diag = diag(n2)*sigma.u.v
Sig = rbind(cbind(Sig.upper.left, Sig.off.diag), cbind(Sig.off.diag, Sig.lower.right))
Sig.corr = cov2cor(Sig)

offdiag.to.vector = function(matrix) {
     n = nrow(matrix)
     vec = rep(0, (n - 1) * n / 2)
     pos = matrix(0, (n-1) * n / 2, 2)
     ind = 1
     for (i in (1:(n-1))) {
       for (j in ((i+1):n)) {
         vec[ind] = matrix[i, j]
         pos[ind, ] = c(i, j)
         ind = ind + 1
       }
     }
     list(vector = vec, position.mat = pos)
   }

On 8/19/09 10:29 AM, "Afshartous, David" <DAfshartous at med.miami.edu> wrote:

Fabian,

RE defining weights, see section 5.2.1 in Pinheiro & Bates (Mixed Effects Models in S and S-Plus).
RE correlation structures, see section 5.3.3.   While these sections are with respect to the univariate mixed model, note that the multivariate mixed model is analogous to the univariate mixed model that is stratified per some grouping such as say gender.  For instance, as far as the data structure goes, there is no difference if we have 2 variables on 10 subjects over time versus  1 variable on 10 subjects over time for males and females.   Of course, the different scenarios will most likely lead to different model assumptions.

While I haven't used the corrrelation = corSymm for a multivariate mixed model in nlme, I have used varIdent to specify different error variances for the different response variables in the multivariate setting.  Just set up your model as specified in the Doran & Lockwood reference below, and include weights = varIdent(form ~ weights = Ind), where Ind is an indicator variable that identifies each separate response variable. What I would like to do is have the error variances also correlated across the different response variables, but only when they are both measured at the same time point.  I'll check this out again and let you know if this and your structure is possible.

Although it's respect to SAS, a good reference that provides insight into the various assumptions in the multivariate mixed model and the resulting structure of the covariance matrices is:
On the use of PROC MIXED to Estimate Correlation in the Presence of Repeated Measures by Hamlett, Ryan, and Wolfinger (comes up via google)

Cheers,
David

On 8/18/09 7:14 AM, "Mollet, Fabian" <Fabian.Mollet at wur.nl> wrote:

Dear nlme expert

We need two pieces of information about the fitting of a nlme model which we cannot extract from the R help files and would be most grateful if you could help us. We fit an energy allocation growth model with 4 parameters to individual growth curves using the nlme routine. We thus have repeated age and size measurements of individuals and therefore allow for random individual effects (i.e. the data is grouped by individual).

1)      Because the sampling of these individuals was size stratified we have to account for the representation of the individual in the true size distribution by statistical weighting. The statistical weight would thus differ across individuals but be the same over the repeated measurements of each individual (to which the random effects apply) and should be somehow multiplied by the residuals of the repeated measurements of each individual. We guess we need to use the varClasses argument but it does not seem clear in the R help files to which level the statistical weights would apply. Could you please tell us how to define the statistical weights on the level of the random effects, i.e. on the level of the individual? varIdent?

2)      We furthermore want to analyze the results of the 4 estimated parameters over time using the lme routine and have thus now 1 row per individual (comprising of the 4 parameters, a time variable and others). Because the 4 parameters are correlated we intend to analyze this multivariate outcome by "flagging" the response by using a dummy coding for the 4 parameters and the time variable as is e.g. described in Doran and Lockwood (2006) p. 223-225 (resulting in 16 rows per individual). Since we want to follow the evolution of the correlation between the 4 parameters over time we would like to make no assumptions on the correlation structure of the errors. We guess we therefore have to use the correlation=corSymm argument. However, the same weighting would apply as in 1) above to the individual and we are therefore not sure again how to define the statistical weights in this case and what this would imply for the error correlation structure. Could you give us a guidance?

Your help is most appreciated and we thank you very much in advance!

Kind regards

Fabian Mollet

Doran, H. C. and Lockwood, J. R. 2006. Fitting value-added models in R. - Journal of Educational and Behavioral Statistics 31: 205-230.

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models