1
$\begingroup$

I am trying to simulate some data from a linear fixed effects model with correlated errors. The model is quite simple, inspired in a pre-post study.

$$Y_i|X_i=x_i\sim N(x_i\beta,\Sigma_i)$$

where $Y_i=(y_{i,pre},y_{i,post})^T$, $x_i=\begin{pmatrix} 1 & 0 \\ 1 & 1 \\ \end{pmatrix}$, $\beta=(\beta_0,\beta_1)^T$ and

$$\Sigma_i=\begin{pmatrix} \sigma_{pre}^2 & \sigma_{pre,post} \\ \sigma_{pre,post} & \sigma_{post}^2 \\ \end{pmatrix}$$

$i$ is indexing the subject from which there are two measures available: pre and post measures.

I am ussing the fact that $\Sigma_i=\Lambda_iR_i\Lambda_i$ where $\Lambda_i=diag(\sigma_{pre},\sigma_{post})$ and $R_i$ is the correlation matrix, which in this case depends only on $\rho_{pre,post}$. Note that $x_i$ and $\Sigma_i$ are equal for all $i$.

Finally a general formulation equation is stated for the model:

$$Y|X\sim N(X\beta,\Sigma)$$

where $Y=(Y_1,\cdots,Y_n)^T$, $X=1_{n \times1} \bigotimes x_i$ and $\Sigma=I_{n \times n}\bigotimes \Sigma_i$ where $\bigotimes$ is the Kronecker Product.

Now let us give some arbitrary values for the model parameters:

$\beta_0=100$, $\beta_1=-30$, $\rho_{pre,post}= -0.8$, and $\sigma_{pre}=\sigma_{post}=9$

Next I present the R code to simulate some data from this model:

library(pacman)
p_load(tidyverse)
p_load(MASS)
p_load(nlme)
p_load(mvtnorm)

##number of subjects
my_n<-500
##correlation between measures inside each subject
my_rho<- -0.8
## sd inside each subject
my_sd<-9
#subject level correlation matrix
sl_mcor<-matrix(c(1,my_rho,my_rho,1),nrow = 2,byrow=T)
#subject level variance (diagonal) matrix
sl_diag_v<-matrix(c(my_sd,0,0,my_sd),nrow = 2,byrow=T)
#subject level covariance matrix
sl_mcov<-sl_diag_v%*%sl_mcor%*%sl_diag_v
#general correlation matrix
mcov<-kronecker(diag(my_n),sl_mcov)
#design matrix for pre-post measures for each subject
z0<-expand.grid(id=(1:my_n) %>% factor,trt=c("pre","post")) %>% as.data.frame()
# fixed effects parameters
my_bts<-c(100,-30)
#model matrix
z0_mat <- model.matrix(~ trt, z0)
# expected values
my_mu<-z0_mat%*%my_bts
#simulate data from Multivariate Normal dist with my_mu and mcov
y<-rmvnorm(1,my_mu,mcov)
z0$y<-y %>% t
#descriptive
z0 %>% ggplot(aes(x=trt,y=y))+geom_boxplot()

#fitting model
m0<-gls(y~trt,correlation=corCompSymm(form=~1|id),data=z0)
summary(m0)
#> Generalized least squares fit by REML
#>   Model: y ~ trt 
#>   Data: z0 
#>        AIC      BIC    logLik
#>   7257.624 7277.247 -3624.812
#> 
#> Correlation Structure: Compound symmetry
#>  Formula: ~1 | id 
#>  Parameter estimate(s):
#>       Rho 
#> 0.1305259 
#> 
#> Coefficients:
#>                 Value Std.Error   t-value p-value
#> (Intercept) 100.02863 0.4081594 245.07248       0
#> trtpost     -29.84501 0.5382365 -55.44962       0
#> 
#>  Correlation: 
#>         (Intr)
#> trtpost -0.659
#> 
#> Standardized residuals:
#>          Min           Q1          Med           Q3          Max 
#> -2.772941799 -0.698851160 -0.003973076  0.710269078  2.923039822 
#> 
#> Residual standard error: 9.126722 
#> Degrees of freedom: 1000 total; 998 residual
intervals(m0,which = "var-cov")
#> Approximate 95% confidence intervals
#> 
#>  Correlation structure:
#>          lower      est.     upper
#> Rho 0.04368246 0.1305259 0.2154116
#> 
#>  Residual standard error:
#>    lower     est.    upper 
#> 8.731751 9.126722 9.539558

Created on 2022-09-19 by the reprex package (v2.0.1)

The concerning issue here is that I am having a hard time recovering the correlation parameter ($\rho_{pre,post}= -0.8$). Estimation for this parameter is quite far from the true parameter set in the simulations. Other parameters are recovered much better than this. What could be the explanation for this? (set aside multiple runs of the simulation and checking the distribution of estimations which I have done already but, this is not shown in the code). The fact that only two measures per subject are available could be a cause for this? is sample size an issue here? $n=500$ is not enough?

As a side question related to this, in the beginning, I was using a positive value for $\rho_{pre,post}= 0.8$ but then I realized that this could be not realistically compatible with the fixed effects parameters $\beta=(100,-30)$, because this fixed effects parameters mostly for sure would induce a negative correlation between $y_{i,pre}$ and $y_{i,post}$ for a given subject. Thus, I am having some doubts about the possible effects on the fixed effects component of the model on the correlation structure itself.

Edit: on the last thought, if we think that residuals are correlated (rather than observations on variable $y$), then fixed effects parameters shouldn't have any effects on the correlation structure. I mean after all $cor(y_{i,pre},y_{i,post})=cor(\epsilon_{i,pre},\epsilon_{i,post})$

$\endgroup$
1
  • $\begingroup$ (+1, generally well-posed question) I fixed two mistakes and one inexact expression in your model formulation. Your $\LaTeX$ code could still benefit from some clean-up. What causes your problem is just a small programming error, not an estimation issue (see my answer). $\endgroup$
    – statmerkur
    Commented Sep 20, 2022 at 9:11

1 Answer 1

2
$\begingroup$

Note that expand.grid() generates a data.frame in which the first factor varies fastest.
Therefore, you need

z0 <- expand.grid(trt = c("pre", "post"), id = 1:my_n %>% factor)

to match the structure of your covariance matrix mcov.

$\endgroup$
3
  • $\begingroup$ indeed the correlation matrix does not fit the structure of the data, I did not care about it too much since I was using the model.matrix function to specify the mean structure of the model, but forgot to take care in the same way for the correlation matrix. Is a good exercise to work this out. ty for your feedback. $\endgroup$ Commented Sep 20, 2022 at 12:48
  • $\begingroup$ @NicolasMolano You're welcome. The way you've specified the statistical model suggests that your code generates a wrong representation of the design matrix $X$, not of the covariance matrix $\Sigma$. That's why I suggested to change the way z0 is generated. $\endgroup$
    – statmerkur
    Commented Sep 20, 2022 at 14:49
  • $\begingroup$ solved, covariance matrix construction for any layout of z0: Z<-model.matrix(~-1+ id, z0);C0<-Z%*%t(Z);mcorr<-(my_rho)*C0+(1-my_rho)*diag(nrow=nrow(C0));mcov<-diag(my_sd,nrow=nrow(C0))%*%mcorr%*%diag(my_sd,nrow=nrow(C0)) just in case any one finds this useful $\endgroup$ Commented Sep 20, 2022 at 20:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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