Correlated Data

Sometimes it is desirable to simulate correlated data from a correlation matrix directly. For example, a simulation might require two random effects (e.g. a random intercept and a random slope). Correlated data like this could be generated using the defData functionality, but it may be more natural to do this with genCorData or addCorData. Currently, simstudy can only generate multivariate normal using these functions.

genCorData requires the user to specify a mean vector mu, a single standard deviation or a vector of standard deviations sigma, and either a correlation matrix corMatrix or a correlation coefficient rho and a correlation structure corsrt. Here are a few examples:

# specifying a specific correlation matrix C
C <- matrix(c(1, 0.7, 0.2, 0.7, 1, 0.8, 0.2, 0.8, 1), nrow = 3)
C
##      [,1] [,2] [,3]
## [1,]  1.0  0.7  0.2
## [2,]  0.7  1.0  0.8
## [3,]  0.2  0.8  1.0
set.seed(282726)

# generate 3 correlated variables with different location and scale for each
# field
dt <- genCorData(1000, mu = c(4, 12, 3), sigma = c(1, 2, 3), corMatrix = C)
dt
## Key: <id>
##          id       V1       V2        V3
##       <int>    <num>    <num>     <num>
##    1:     1 4.125728 12.92567  3.328106
##    2:     2 4.712100 14.26502  8.876664
##    3:     3 4.990881 14.44321  5.322747
##    4:     4 4.784358 14.86861  8.129774
##    5:     5 4.930617 11.11235 -1.400923
##   ---                                  
##  996:   996 2.983723 13.61509  8.773969
##  997:   997 2.852707 10.43317  3.811047
##  998:   998 3.856643 13.17697  4.720628
##  999:   999 4.738479 12.64438  2.979415
## 1000:  1000 5.766867 13.51827  1.693172
# estimate correlation matrix
dt[, round(cor(cbind(V1, V2, V3)), 1)]
##     V1  V2  V3
## V1 1.0 0.7 0.2
## V2 0.7 1.0 0.8
## V3 0.2 0.8 1.0
# estimate standard deviation
dt[, round(sqrt(diag(var(cbind(V1, V2, V3)))), 1)]
##  V1  V2  V3 
## 0.9 1.9 3.0
# generate 3 correlated variables with different location but same standard
# deviation and compound symmetry (cs) correlation matrix with correlation
# coefficient = 0.4.  Other correlation matrix structures are 'independent'
# ('ind') and 'auto-regressive' ('ar1').

dt <- genCorData(1000, mu = c(4, 12, 3), sigma = 3, rho = 0.4, corstr = "cs", cnames = c("x0",
    "x1", "x2"))
dt
## Key: <id>
##          id        x0        x1         x2
##       <int>     <num>     <num>      <num>
##    1:     1 7.1160161 14.294748  4.0251237
##    2:     2 3.5429823  8.299333  4.5620657
##    3:     3 2.5590428 10.660403  2.5805860
##    4:     4 5.9808506 12.457614  2.0287775
##    5:     5 3.7210289 15.003835  7.4425421
##   ---                                     
##  996:   996 0.3996175  8.104629  5.5241810
##  997:   997 1.4299019 11.311426 -0.6144622
##  998:   998 3.3079075 11.909745 -0.7375013
##  999:   999 3.7934154 10.515881  2.6021325
## 1000:  1000 5.6413141 13.513672  7.5321371
# estimate correlation matrix
dt[, round(cor(cbind(x0, x1, x2)), 1)]
##     x0  x1  x2
## x0 1.0 0.5 0.4
## x1 0.5 1.0 0.4
## x2 0.4 0.4 1.0
# estimate standard deviation
dt[, round(sqrt(diag(var(cbind(x0, x1, x2)))), 1)]
##  x0  x1  x2 
## 2.9 3.0 3.0

The new data generated by genCorData can be merged with an existing data set. Alternatively, addCorData will do this directly:

# define and generate the original data set
def <- defData(varname = "x", dist = "normal", formula = 0, variance = 1, id = "cid")
dt <- genData(1000, def)

# add new correlate fields a0 and a1 to 'dt'
dt <- addCorData(dt, idname = "cid", mu = c(0, 0), sigma = c(2, 0.2), rho = -0.2,
    corstr = "cs", cnames = c("a0", "a1"))

dt
## Key: <cid>
##         cid          x          a0          a1
##       <int>      <num>       <num>       <num>
##    1:     1 -0.4707940  0.97711194 -0.09127123
##    2:     2 -1.8723668  2.70498417 -0.27102780
##    3:     3  1.3347964 -5.15138578 -0.12289563
##    4:     4 -0.1685203  1.04733271 -0.04400129
##    5:     5 -1.6308055  0.39516494  0.06640973
##   ---                                         
##  996:   996 -0.5244473 -0.61310062  0.27520456
##  997:   997  1.4965903  1.32673834  0.47458481
##  998:   998  0.1015744 -0.09821567  0.06440723
##  999:   999 -0.5788317  0.16870967 -0.03890117
## 1000:  1000 -1.6175613  3.61182553 -0.46220263
# estimate correlation matrix
dt[, round(cor(cbind(a0, a1)), 1)]
##      a0   a1
## a0  1.0 -0.2
## a1 -0.2  1.0
# estimate standard deviation
dt[, round(sqrt(diag(var(cbind(a0, a1)))), 1)]
##  a0  a1 
## 2.0 0.2

Correlated data: additional distributions

Two additional functions facilitate the generation of correlated data from binomial, poisson, gamma, and uniform distributions: genCorGen and addCorGen.

genCorGen is an extension of genCorData. These functions draw on copula-based methods to generate the data. (This Wikipedia page provides a general introduction and copula-based modeling can be conducted in R using package copula.) In the first example, we are generating data from a multivariate Poisson distribution. We start by specifying the mean of the Poisson distribution for each new variable, and then we specify the correlation structure, just as we did with the normal distribution.

l <- c(8, 10, 12) # lambda for each new variable

dx <- genCorGen(1000, nvars = 3, params1 = l, dist = "poisson", rho = .3, corstr = "cs", wide = TRUE)
dx
## Key: <id>
##          id    V1    V2    V3
##       <int> <num> <num> <num>
##    1:     1     5    16    13
##    2:     2     9     9     6
##    3:     3     7    11    18
##    4:     4    11    14    12
##    5:     5    10     8    15
##   ---                        
##  996:   996     3     2     5
##  997:   997     6    14    11
##  998:   998     6     8    12
##  999:   999    10    12    11
## 1000:  1000     9     9    12
round(cor(as.matrix(dx[, .(V1, V2, V3)])), 2)
##     V1   V2   V3
## V1 1.0 0.30 0.30
## V2 0.3 1.00 0.24
## V3 0.3 0.24 1.00

We can also generate correlated binary data by specifying the probabilities:

genCorGen(1000, nvars = 3, params1 = c(.3, .5, .7), dist = "binary", rho = .8, corstr = "cs", wide = TRUE)
## Key: <id>
##          id    V1    V2    V3
##       <int> <num> <num> <num>
##    1:     1     0     1     1
##    2:     2     0     1     1
##    3:     3     0     0     1
##    4:     4     1     1     1
##    5:     5     0     1     1
##   ---                        
##  996:   996     0     0     1
##  997:   997     1     1     1
##  998:   998     0     0     0
##  999:   999     0     0     1
## 1000:  1000     0     0     0

The gamma distribution requires two parameters - the mean and dispersion. (These are converted into shape and rate parameters more commonly used.)

dx <- genCorGen(1000, nvars = 3, params1 = l, params2 = c(1,1,1), dist = "gamma", rho = .7, corstr = "cs", wide = TRUE, cnames="a, b, c")
dx
## Key: <id>
##          id            a          b           c
##       <int>        <num>      <num>       <num>
##    1:     1 4.137889e+00  1.9736693  5.73317661
##    2:     2 6.230611e-04  0.1790216  0.01098133
##    3:     3 9.554613e+00 21.3956071 30.07914569
##    4:     4 1.053229e+01  6.8598915  7.47104860
##    5:     5 2.556925e+01 22.8862611 17.32239223
##   ---                                          
##  996:   996 2.635737e+00  0.9269903  1.22333746
##  997:   997 1.638308e+00 12.0692638 13.01943662
##  998:   998 3.492819e+00  4.1504352  2.37403911
##  999:   999 9.336809e+00 21.2184483 25.17933311
## 1000:  1000 2.044966e+01 32.3326247 23.81715119
round(cor(as.matrix(dx[, .(a, b, c)])), 2)
##      a    b    c
## a 1.00 0.63 0.67
## b 0.63 1.00 0.67
## c 0.67 0.67 1.00

These data sets can be generated in either wide or long form. So far, we have generated wide form data, where there is one row per unique id. Now, we will generate data using the long form, where the correlated data are on different rows, so that there are repeated measurements for each id. An id will have multiple records (i.e. one id will appear on multiple rows):

dx <- genCorGen(1000, nvars = 3, params1 = l, params2 = c(1,1,1), dist = "gamma", rho = .7, corstr = "cs", wide = FALSE, cnames="NewCol")
dx
## Key: <id>
##          id period     NewCol
##       <int>  <num>      <num>
##    1:     1      0 0.08868527
##    2:     1      1 0.17558015
##    3:     1      2 0.35553817
##    4:     2      0 2.41522425
##    5:     2      1 0.99489378
##   ---                        
## 2996:   999      1 4.62541703
## 2997:   999      2 0.73199287
## 2998:  1000      0 3.52197152
## 2999:  1000      1 7.43262675
## 3000:  1000      2 8.36619208

addCorGen allows us to create correlated data from an existing data set, as one can already do using addCorData. In the case of addCorGen, the parameter(s) used to define the distribution are created as a field (or fields) in the dataset. The correlated data are added to the existing data set. In the example below, we are going to generate three sets (poisson, binary, and gamma) of correlated data with means that are a function of the variable xbase, which varies by id.

First we define the data and generate a data set:

def <- defData(varname = "xbase", formula = 5, variance = .2, dist = "gamma", id = "cid")
def <- defData(def, varname = "lambda", formula = ".5 + .1*xbase", dist="nonrandom", link = "log")
def <- defData(def, varname = "p", formula = "-2 + .3*xbase", dist="nonrandom", link = "logit")
def <- defData(def, varname = "gammaMu", formula = ".5 + .2*xbase", dist="nonrandom", link = "log")
def <- defData(def, varname = "gammaDis", formula = 1, dist="nonrandom")

dt <- genData(10000, def)
dt
## Key: <cid>
##          cid    xbase   lambda         p  gammaMu gammaDis
##        <int>    <num>    <num>     <num>    <num>    <num>
##     1:     1 1.546326 1.924435 0.1771026 2.246257        1
##     2:     2 5.689908 2.912439 0.4272628 5.144775        1
##     3:     3 5.059867 2.734604 0.3817705 4.535672        1
##     4:     4 4.599528 2.611573 0.3497493 4.136730        1
##     5:     5 2.402442 2.096447 0.2176749 2.665758        1
##    ---                                                    
##  9996:  9996 3.610769 2.365707 0.2856166 3.394491        1
##  9997:  9997 4.984305 2.714019 0.3764348 4.467643        1
##  9998:  9998 5.122724 2.751847 0.3862310 4.593052        1
##  9999:  9999 3.393940 2.314964 0.2725312 3.250432        1
## 10000: 10000 7.722561 3.568895 0.5785365 7.725390        1

The Poisson distribution has a single parameter, lambda:

dtX1 <- addCorGen(dtOld = dt, idvar = "cid", nvars = 3, rho = .1, corstr = "cs",
                    dist = "poisson", param1 = "lambda", cnames = "a, b, c")
dtX1
## Key: <cid>
##          cid    xbase   lambda         p  gammaMu gammaDis     a     b     c
##        <int>    <num>    <num>     <num>    <num>    <num> <num> <num> <num>
##     1:     1 1.546326 1.924435 0.1771026 2.246257        1     2     2     2
##     2:     2 5.689908 2.912439 0.4272628 5.144775        1     1     0     2
##     3:     3 5.059867 2.734604 0.3817705 4.535672        1     4     0     2
##     4:     4 4.599528 2.611573 0.3497493 4.136730        1     3     1     3
##     5:     5 2.402442 2.096447 0.2176749 2.665758        1     5     4     1
##    ---                                                                      
##  9996:  9996 3.610769 2.365707 0.2856166 3.394491        1     1     2     2
##  9997:  9997 4.984305 2.714019 0.3764348 4.467643        1     3     2     4
##  9998:  9998 5.122724 2.751847 0.3862310 4.593052        1     7     2     1
##  9999:  9999 3.393940 2.314964 0.2725312 3.250432        1     6     0     4
## 10000: 10000 7.722561 3.568895 0.5785365 7.725390        1    12     1     3

The Bernoulli (binary) distribution has a single parameter, p:

dtX2 <- addCorGen(dtOld = dt, idvar = "cid", nvars = 4, rho = .4, corstr = "ar1",
                    dist = "binary", param1 = "p")
dtX2
## Key: <cid>
##          cid    xbase   lambda         p  gammaMu gammaDis    V1    V2    V3
##        <int>    <num>    <num>     <num>    <num>    <num> <num> <num> <num>
##     1:     1 1.546326 1.924435 0.1771026 2.246257        1     0     0     1
##     2:     2 5.689908 2.912439 0.4272628 5.144775        1     0     0     0
##     3:     3 5.059867 2.734604 0.3817705 4.535672        1     1     0     1
##     4:     4 4.599528 2.611573 0.3497493 4.136730        1     0     0     0
##     5:     5 2.402442 2.096447 0.2176749 2.665758        1     0     1     0
##    ---                                                                      
##  9996:  9996 3.610769 2.365707 0.2856166 3.394491        1     0     0     0
##  9997:  9997 4.984305 2.714019 0.3764348 4.467643        1     0     1     1
##  9998:  9998 5.122724 2.751847 0.3862310 4.593052        1     0     0     0
##  9999:  9999 3.393940 2.314964 0.2725312 3.250432        1     0     0     0
## 10000: 10000 7.722561 3.568895 0.5785365 7.725390        1     0     0     1
##           V4
##        <num>
##     1:     0
##     2:     1
##     3:     1
##     4:     1
##     5:     0
##    ---      
##  9996:     0
##  9997:     0
##  9998:     1
##  9999:     0
## 10000:     0

The Gamma distribution has two parameters - in simstudy the mean and dispersion are specified:

dtX3 <- addCorGen(dtOld = dt, idvar = "cid", nvars = 4, rho = .4, corstr = "cs",
                  dist = "gamma", param1 = "gammaMu", param2 = "gammaDis")
dtX3
## Key: <cid>
##          cid    xbase   lambda         p  gammaMu gammaDis          V1
##        <int>    <num>    <num>     <num>    <num>    <num>       <num>
##     1:     1 1.546326 1.924435 0.1771026 2.246257        1 4.231680194
##     2:     2 5.689908 2.912439 0.4272628 5.144775        1 6.787710358
##     3:     3 5.059867 2.734604 0.3817705 4.535672        1 3.707544002
##     4:     4 4.599528 2.611573 0.3497493 4.136730        1 0.348766326
##     5:     5 2.402442 2.096447 0.2176749 2.665758        1 2.493923583
##    ---                                                                
##  9996:  9996 3.610769 2.365707 0.2856166 3.394491        1 0.002486867
##  9997:  9997 4.984305 2.714019 0.3764348 4.467643        1 3.210241942
##  9998:  9998 5.122724 2.751847 0.3862310 4.593052        1 9.556894110
##  9999:  9999 3.393940 2.314964 0.2725312 3.250432        1 1.349413306
## 10000: 10000 7.722561 3.568895 0.5785365 7.725390        1 2.404109193
##               V2         V3         V4
##            <num>      <num>      <num>
##     1: 1.2176380  0.4222348  2.1937488
##     2: 4.6349836  6.5642077  3.1765033
##     3: 1.5053746  1.1395938  1.0282041
##     4: 7.7065647 15.1843651  1.5413112
##     5: 0.6724540  2.0454487  1.5117235
##    ---                                
##  9996: 0.4825269  0.6716621  0.6721581
##  9997: 0.2103485  4.9869686  3.8444291
##  9998: 3.5482244 10.5174231  9.9577729
##  9999: 2.3406277  1.7019004  0.8970987
## 10000: 2.7987448  3.6667012 12.1680528

If we have data in long form (e.g. longitudinal data), the function will recognize the structure:

def <- defData(varname = "xbase", formula = 5, variance = .4, dist = "gamma", id = "cid")
def <- defData(def, "nperiods", formula = 3, dist = "noZeroPoisson")

def2 <- defDataAdd(varname = "lambda", formula = ".5+.5*period + .1*xbase", dist="nonrandom", link = "log")

dt <- genData(1000, def)

dtLong <- addPeriods(dt, idvars = "cid", nPeriods = 3)
dtLong <- addColumns(def2, dtLong)

dtLong
## Key: <timeID>
##         cid period    xbase nperiods timeID    lambda
##       <int>  <int>    <num>    <num>  <int>     <num>
##    1:     1      0 7.053471        2      1  3.337917
##    2:     1      1 7.053471        2      2  5.503295
##    3:     1      2 7.053471        2      3  9.073400
##    4:     2      0 2.185136        3      4  2.051382
##    5:     2      1 2.185136        3      5  3.382157
##   ---                                                
## 2996:   999      1 9.702454        1   2996  7.172436
## 2997:   999      2 9.702454        1   2997 11.825348
## 2998:  1000      0 3.044209        4   2998  2.235402
## 2999:  1000      1 3.044209        4   2999  3.685554
## 3000:  1000      2 3.044209        4   3000  6.076452
### Generate the data 

dtX3 <- addCorGen(dtOld = dtLong, idvar = "cid", nvars = 3, rho = .6, corstr = "cs",
                  dist = "poisson", param1 = "lambda", cnames = "NewPois")
dtX3
## Key: <cid>
##         cid period    xbase nperiods timeID    lambda NewPois
##       <int>  <int>    <num>    <num>  <int>     <num>   <num>
##    1:     1      0 7.053471        2      1  3.337917       5
##    2:     1      1 7.053471        2      2  5.503295       6
##    3:     1      2 7.053471        2      3  9.073400      11
##    4:     2      0 2.185136        3      4  2.051382       2
##    5:     2      1 2.185136        3      5  3.382157       1
##   ---                                                        
## 2996:   999      1 9.702454        1   2996  7.172436      10
## 2997:   999      2 9.702454        1   2997 11.825348      11
## 2998:  1000      0 3.044209        4   2998  2.235402       4
## 2999:  1000      1 3.044209        4   2999  3.685554       7
## 3000:  1000      2 3.044209        4   3000  6.076452       9

We can fit a generalized estimating equation (GEE) model and examine the coefficients and the working correlation matrix. They match closely to the data generating parameters:

geefit <- gee(NewPois ~ period + xbase, data = dtX3, id = cid, family = poisson, corstr = "exchangeable")
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
## (Intercept)      period       xbase 
##   0.4915131   0.4929285   0.1030995
round(summary(geefit)$working.correlation, 2)
##      [,1] [,2] [,3]
## [1,] 1.00 0.56 0.56
## [2,] 0.56 1.00 0.56
## [3,] 0.56 0.56 1.00