I wish to conduct a path analysis (graph below) on data which are repeated for some individuals (but not all individuals, nor all at the same time). Is it possible to model variance structure with for example an individual variance or temporal variance (between years) in a path analysis? If it is, in which software and how can I code the different sources of error variances? I am quite familiar with the SEM package in R, so a method to model longitudinal data in SEM would be ideal.

Bonus question: Is it possible to model different variance structures for the different variables in the model (as is possible for multivariate regression)?

Thanks for your help!

Path analysis graph


The solution using the mcglm package.

Extra packages


Simulating the data set



env1 <- rnorm(100)
env2 <- rnorm(100)

Trait 1

par1 <- 5
par2 <- 2
var1 <- 0.5
mean1 <- par1*env1 + par2*env2

Trait 2

par3 = -1
par4 = 3
var2 = 1
mean2 <- par3*env1 + par4*env2

Covariance between trait1 and trait2

Sigma <- matrix(c(1,0.8,0.8,2), 2, 2)
trait12 <- matrix(NA, ncol = 2, nrow = 100)
for(i in 1:100) {
trait12[i,] <- as.numeric(rmvnorm(n = 1, 
mean = c(mean1[i], mean2[i]),sigma = Sigma))}
trait1 <- trait12[,1]
trait2 <- trait12[,2]

Trait 3

par5 = 2
par6 = -0.3
par7 = 0.10
par8 = 0.5

Note that in the case of binary trait, this variance is 1. var3 = 0.25

Binary trait

#mean3 <- par5*env1 + par6*env2 + par7*trait1 + par8*trait2
#trait3 <- rnorm(100, mean = mean3, sd = sqrt(var3))
mean3 <- mc_link_function(beta = c(par5, par6, par7, par8), 
X = model.matrix(~ env1 + env2 + trait1 + trait2 -1),
link = "logit", offset = NULL)$mu
trait3 <- rbinom(100, size = 1, prob = mean3)


par9 = 1.5
par10 = -0.5
par11 = 3.2
var4 = 0.8
mean.out <- par9*trait1 + par10*trait2 + par11*trait3
outcome <- rnorm(100, mean.out, sd = sqrt(var4))

Your data should be like this data <- data.frame(outcome, trait1, trait2, trait3, env1, env2)

Linear predictor

form.out <- outcome ~ trait1 + trait2 + trait3 -1
form.t1 <- trait1 ~ env1 + env2 -1
form.t2 <- trait2 ~ env1 + env2 -1
form.t3 <- trait3 ~ env1 + env2 + trait1 + trait2 -1

Matrix linear predictor for your data just a diagonal matrix

Z0 <- mc_id(data)

To deal with longitudinal or repeated measures, we have to include a non-diagonal matrix in the matrix linear predictor. Suppose, we have 10 subjects for each one 10 replications. First, I will create the index

data$id <- rep(1:10, each = 10)

I will use a compound symmetry structure very popular in repeated measures basically, it is a matrix of ones

Z1.list <- list()
for(i in 1:10) {
ONE <- rep(1,10)
Z1.list[[i]] <- tcrossprod(ONE)}
Z1 <- bdiag(Z1.list)

Basically, this matrix describes your repeated measures structure. Now, just include it in the matrix linear predictor. Note that, I did not simulate the data using it, so the parameter associated with this structure should be close to zero. I have some functions that help to construct this matrix, but in your case, I think you have some irregular and different levels, so it is better you construct your matrices. Just use the Matrix package to do it!!!

We can create the same matrix using

Z12 <- mc_mixed(~ 0 + as.factor(id), data = data)

but if your structure is not balance, this function are now very usefull!

Fitting the model


Fit without repeated measures structure

fit1 <- mcglm(linear_pred = c(form.out, form.t1, form.t2,form.t3),
matrix_pred = list(c(Z0), c(Z0), c(Z0), c(Z0)), 
link = c("identity", "identity", "identity", "logit"),
variance = c("constant", "constant", "constant", "binomialP"),
data = data, 
control_algorithm = list(tol = 1e-03, max_iter = 100, 
tunning = 1, verbose = TRUE))
plot(fit1, type = "algorithm") # OK

In general, if you include many components in the matrix linear predictor the model is very hard to fit. But, anyway to include extra components in the matrix linear predictor is quite easy see below.

Fit with repeated measures structure

fit2 <- mcglm(linear_pred = c(form.out, form.t1, form.t2, form.t3), 
matrix_pred = list(c(Z0, Z1), c(Z0, Z1), c(Z0, Z1), c(Z0, Z1)), 
link = c("identity", "identity", "identity", "logit"),
variance = c("constant", "constant", "constant", "binomialP"),
data = data, 
control_algorithm = list(tol = 1e-03, max_iter = 100, 
tunning = 1, verbose = TRUE))
plot(fit2, type = "algorithm") # OK
  • $\begingroup$ I managed to run the model with my data, but I still get the same message error: "Cholesky factorization failed". I expect it cannot converge due to my data structure, but I do not have diagnostic tools to understand whether more iterations could help. $\endgroup$ – Charlotte R Sep 21 '16 at 7:56
  • $\begingroup$ Moreover, the argument verbose = TRUE should print covariance parameters for each iteration, but here is a sample of what I get at last iteration: [1] -0.5032 0.4572 0.1339 -0.2076 0.2316 -1.1555 0.8595 -1.3032 1.0275 [10] -0.7348 19.7512 -4.0867 -0.0286 0.4096 [1] 5128.7993 11495.2547 -6209.2595 7510.7570 -4038.8713 -9056.8767 [7] -582.7147 -894.6413 109.9523 -191.7047 3217.2540 291.2220 [13] -275.5062 1643.8223 $\endgroup$ – Charlotte R Sep 21 '16 at 7:57
  • $\begingroup$ And with no binomial variable, but an individual random effects, I get again "Cholesky factorization failed" $\endgroup$ – Charlotte R Sep 21 '16 at 8:06
  • $\begingroup$ With only the diagonal matrix (no random effects), the model does not converge either: Error in solve.default(cov_temp$Sensitivity, cov_temp$Score) : system is computationally singular: reciprocal condition number = 2.47436e-48 However, it does not seem to be a problem with the design matrix itself, since the model does converge OK when fitting a gaussian variance for the binomial variable (without random effects). $\endgroup$ – Charlotte R Sep 21 '16 at 8:31

Based only on this error is difficult to figure out what is going on. This is the unique error possible to occur when using mcglm. It is just telling you that the variance-covariance matrix is not positive definite. Assuming, that your model is correct specified and you data support your model, if you control the step-length through the tunning argument it should converge. I suggest to use tunning = 0.5 or tunning = 0.1.


  • $\begingroup$ Reducing the tunning argument to 0.5 got the model converging with a binomial variable but not with the individual variance, whereas a value of 0.1 did the trick. $\endgroup$ – Charlotte R Sep 27 '16 at 11:48
  • $\begingroup$ However, I am not sure this model is suitable for a path analysis: 1) Estimating simultaneously the regression of Y on X and the correlation between X and Y produces spurious results, e.g. a strong positive correlation but a negative effect, when Y is actually computed as a mathematical linear combination of X and other variables; 2) The covariance structure does not reflect the variance structure, which seems odd to me, i.e. only residual covariance can be estimated, even if a nested variance structure is implemented. $\endgroup$ – Charlotte R Sep 29 '16 at 9:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.