Problem with hlme::predictY

[daavic]

Thanks very much for your lcmm package. I find it very useful for building GMM models with hlme()

I can't understand what predictY is doing after I build the GMM models.

Using time <- seq(from = min(dat.use$age), to = max(dat.use$age), length.out = 100)

time ranges from 7 to 50

I convert this to splines (the model was built using splines as the outcome is highly non-linear) spline.pred <- bs(age.eps, degree = 4, df = 6)

Then I give the splines to predictY() and I get my one hundred predictions

I adjust the time variable, time.eps <- time + eps

eps is small ~.4, (I am wanting to calculate derivatives)

I convert this to splines, spline.eps <- bs(time.eps, degree = 4, df = 6)

When I feed this to predictY(), I get exactly the same results as for the original splines without the eps.

I don't understand why this would be the case?

However, if I do time <- seq(from = min(dat.use$age), to = max(dat.use$age), length.out = 200)

The new 100 age inputs, between the original 100, get new predictions.

Am I doing something wrong?

[daavic]

While I think about it, is it possible to get confidence limits for the predictions?

[daavic]

In the covariance matrix for a hlme model obtained with vcov(hlme.model), I notice some of the column names are repeated with a .1 suffix. For example, I have one column with "intercept class1" and a second with "intercept class1.1". Why are they named like this - what do they represent - and which one would I use to formulate a confidence interval for a prediction?

[VivianePhilipps]

Hi,

Here are my answers :

1. with splines, you need to specify the knots, not the degree. See for example #67.
2. the CI for the predictions are obtained with draws=TRUE in predictY.
3. the first intercept refers to the parameter in the latent class membership model, whereas the second intercept refers to the parameter of the mixed model. If you want to check it, the order in vcov(hlme.model) is the same as in hlme.model$best, and the values in hlme.model$best are those printed in summary(hlme.model), where the names are more explicit.

Best,

Viviane

[daavic]

Thanks very much Viviane - much appreciated.

I missed the draws argument in the CRAN document.

Before I got your reply I worked out another way of getting the CI.

get the covariance matrix

var_beta <- vcov(hlme_mod)

subset the rows and cols required

var_beta <- (var_beta)[c(3:23), c(3:23)]

just need to get class2 variances

class2_names <- grep("class2$", rownames(var_beta), value = TRUE)

var_beta.class2 <- var_beta_fixed[class2_names, class2_names]

model matrix s1 etc are covariate names

X <- model.matrix(~ s1+s2+s3+s4+s5+s6, data = newdata )

Variance of predictions

var_y_hat.class2 <- X %% var_beta.class2 %% t(X)

Variance of predictions, standard errors = sqrt

var.class2 <- (diag(var_y_hat.class2))

The draws argument is much better!

Thanks again Viviane for taking the time - lcmm is a very useful package.