## Overview

The Stochastic Process Model (SPM) was developed several decades ago (Woodbury and Manton 1977, A. I. Yashin, Arbeev, Akushevich, et al. (2007)), and applied for analyses of clinical, demographic, epidemiologic longitudinal data as well as in many other studies that relate stochastic dynamics of repeated measures to the probability of end-points (outcomes). SPM links the dynamic of stochastical variables with a hazard rate as a quadratic function of the state variables (A. I. Yashin, Arbeev, Akushevich, et al. 2007). The R-package, “stpm”, is a set of utilities to estimate parameters of stochastic process and modeling survival trajectories and time-to-event outcomes observed from longitudinal studies. It is a general framework for studying and modeling survival (censored) traits depending on random trajectories (stochastic paths) of variables.

## Installation

### Stable version from CRAN

install.packages("stpm")

require(devtools)
devtools::install_github("izhbannikov/stpm")

## Data description

Data represents a typical longitudinal data in form of two datasets: longitudinal dataset (follow-up studies), in which one record represents a single observation, and vital (survival) statistics, where one record represents all information about the subject. Longitudinal dataset cat contain a subject ID (identification number), status (event(1)/censored(0)), time and measurements across the variables.

Below there is an example of clinical data that can be used in stpm and we will discuss the fields later.

Longitudinal table:

##   ID IndicatorDeath Age      DBP      BMI
## 1  1              0  30 80.00000 25.00000
## 2  1              0  32 80.51659 26.61245
## 3  1              0  34 77.78412 29.16790
## 4  1              0  36 77.86665 32.40359
## 5  1              0  38 96.55673 31.92014
## 6  1              0  40 94.48616 32.89139

#### Description of data fields

• ID - subject unique identificatin number.
• IndicatorDeath - 0/1, indicates death of a subject.
• Age - current age of subject at observation.
• DBP, BMI - covariates, here “DBP” represents a diastolic blood pressure, “BMI” a body-mass index.

#### “Short” and “Long” longitudinal data formats

The packate accepts longitudinal data in two formats: “short” and “long”.

##### “Short” format
##   id xi  t       y1
## 1  1  0 30 72.51327
## 2  1  0 31 72.72882
## 3  1  0 32 82.53504
## 4  1  0 33 86.18644
## 5  1  0 34 83.16840
## 6  1  0 35 88.83286
##### “Long” format
##   id xi t1 t2       y1  y1.next
## 1  1  0 30 31 78.26699 73.43093
## 2  1  0 31 32 73.43093 70.44998
## 3  1  0 32 33 70.44998 64.03091
## 4  1  0 33 34 64.03091 62.79530
## 5  1  0 34 35 62.79530 66.77295
## 6  1  0 35 36 66.77295 67.64856

## Discrete- and continuous-time models

There are two main SPM types in the package: discrete-time model (Akushevich, Kulminski, and Manton 2005) and continuous-time model (A. I. Yashin, Arbeev, Akushevich, et al. 2007). Discrete model assumes equal intervals between follow-up observations. The example of discrete dataset is given below.

library(stpm)
data <- simdata_discr(N=10) # simulate data for 10 individuals, "long" format (default)
head(data)
##   id xi t1 t2       y1  y1.next
## 1  1  0 30 31 85.99085 88.61096
## 2  1  0 31 32 88.61096 86.76278
## 3  1  0 32 33 86.76278 93.51249
## 4  1  0 33 34 93.51249 81.05101
## 5  1  0 34 35 81.05101 81.78848
## 6  1  0 35 36 81.78848 80.93441

In this case there are equal intervals between $$t_1$$ and $$t_2$$.

In the continuous-time SPM, in which intervals between observations are not equal (arbitrary or random). The example of such dataset is shown below:

library(stpm)
data <- simdata_cont(N=5, format="short") # simulate data for 5 individuals, "short" format
head(data)
##   id xi        t       y1
## 1  0  0 38.33163 81.38445
## 2  0  0 39.52192 78.85477
## 3  0  0 40.99396 84.38290
## 4  0  0 42.25867 76.14786
## 5  0  0 43.56306 84.89407
## 6  0  0 44.69685 76.07064

### Discrete-time model

The discrete model assumes fixed time intervals between consecutive observations. In this model, $$\mathbf{Y}(t)$$ (a $$k \times 1$$ matrix of the values of covariates, where $$k$$ is the number of considered covariates) and $$\mu(t, \mathbf{Y}(t))$$ (the hazard rate) have the following form:

$$\mathbf{Y}(t+1) = \mathbf{u} + \mathbf{R} \mathbf{Y}(t) + \mathbf{\epsilon}$$

$$\mu (t, \mathbf{Y}(t)) = [\mu_0 + \mathbf{b} \mathbf{Y}(t) + \mathbf{Y}(t)^* \mathbf{Q} \mathbf{Y}(t)] e^{\theta t}$$

Coefficients $$\mathbf{u}$$ (a $$k \times 1$$ matrix, where $$k$$ is a number of covariates), $$\mathbf{R}$$ (a $$k \times k$$ matrix), $$\mu_0$$, $$\mathbf{b}$$ (a $$1 \times k$$ matrix), $$\mathbf{Q}$$ (a $$k \times k$$ matrix) are assumed to be constant in the particular implementation of this model in the R-package stpm. $$\mathbf{\epsilon}$$ are normally-distributed random residuals, $$k \times 1$$ matrix. A symbol ’*’ denotes transpose operation. $$\theta$$ is a parameter to be estimated along with other parameters ($$\mathbf{u}$$, $$\mathbf{R}$$, $$\mathbf{\mu_0}$$, $$\mathbf{b}$$, $$\mathbf{Q}$$).

#### Example

library(stpm)
#Data simulation (200 individuals)
data <- simdata_discr(N=100)
#Estimation of parameters
pars <- spm_discrete(data)
pars
## $dmodel ##$dmodel$theta ## [1] 0.088 ## ##$dmodel$mu0 ## [1] 8.490158093e-05 ## ##$dmodel$b ## [1] -2.062971766e-06 ## ##$dmodel$Q ## [,1] ## [1,] 1.315223336e-08 ## ##$dmodel$u ## [1] 4.289968822 ## ##$dmodel$u.std.err ## (Intercept) ## 0.3539033373 ## ##$dmodel$R ## [,1] ## [1,] 0.9473910422 ## ##$dmodel$R.std.err ## y1_1 ## [1,] 0.004260951529 ## ##$dmodel$Sigma ## [1] 4.992027889 ## ## ##$cmodel
## $cmodel$a
##                [,1]
## [1,] -0.05260895781
##
## $cmodel$f1
##             [,1]
## [1,] 81.54445556
##
## $cmodel$Q
##                 [,1]
## [1,] 1.315223336e-08
##
## $cmodel$f
##             [,1]
## [1,] 78.42667134
##
## $cmodel$b
##             [,1]
## [1,] 4.992027889
##
## $cmodel$mu0
##                 [,1]
## [1,] 4.005576577e-06
##
## $cmodel$theta
## [1] 0.088
##
##
## attr(,"class")
## [1] "spm.discrete"

### Continuous-time model

In the specification of the SPM described in 2007 paper by Yashin and collegaues (A. I. Yashin, Arbeev, Akushevich, et al. 2007) the stochastic differential equation describing the age dynamics of a covariate is:

$$d\mathbf{Y}(t)= \mathbf{a}(t)(\mathbf{Y}(t) -\mathbf{f}_1(t))dt + \mathbf{b}(t)d\mathbf{W}(t), \mathbf{Y}(t=t_0)$$

In this equation, $$\mathbf{Y}(t)$$ (a $$k \times 1$$ matrix) is the value of a particular covariate at a time (age) $$t$$. $$\mathbf{f}_1(t)$$ (a $$k \times 1$$ matrix) corresponds to the long-term mean value of the stochastic process $$\mathbf{Y}(t)$$, which describes a trajectory of individual covariate influenced by different factors represented by a random Wiener process $$\mathbf{W}(t)$$. Coefficient $$\mathbf{a}(t)$$ (a $$k \times k$$ matrix) is a negative feedback coefficient, which characterizes the rate at which the process reverts to its mean. In the area of research on aging, $$\mathbf{f}_1(t)$$ represents the mean allostatic trajectory and $$\mathbf{a}(t)$$ represents the adaptive capacity of the organism. Coefficient $$\mathbf{b}(t)$$ (a $$k \times 1$$ matrix) characterizes a strength of the random disturbances from Wiener process $$\mathbf{W}(t)$$.

The following function $$\mu(t, \mathbf{Y}(t))$$ represents a hazard rate:

$$\mu(t, \mathbf{Y}(t)) = \mu_0(t) + (\mathbf{Y}(t) - \mathbf{f}(t))^* \mathbf{Q}(t) (\mathbf{Y}(t) - \mathbf{f}(t))$$

here $$\mu_0(t)$$ is the baseline hazard, which represents a risk when $$\mathbf{Y}(t)$$ follows its optimal trajectory; $$\mathbf{f}(t)$$ (a $$k \times 1$$ matrix) represents the optimal trajectory that minimizes the risk and $$\mathbf{Q}(t)$$ ($$k \times k$$ matrix) represents a sensitivity of risk function to deviation from the norm.

#### Example

library(stpm)
#Simulate some data for 50 individuals
data <- simdata_cont(N=50)
head(data)
##   id xi          t1          t2          y1     y1.next
## 1  0  0 37.71271511 39.66801034 80.96234637 78.04555132
## 2  0  0 39.66801034 41.27990719 78.04555132 72.88288193
## 3  0  0 41.27990719 42.58912587 72.88288193 70.90928841
## 4  0  0 42.58912587 43.75811857 70.90928841 73.42800255
## 5  0  0 43.75811857 45.74553118 73.42800255 67.82994172
## 6  0  0 45.74553118 47.60762597 67.82994172 75.52968845
#Estimate parameters
# a=-0.05, f1=80, Q=2e-8, f=80, b=5, mu0=2e-5, theta=0.08 are starting values for estimation procedure
pars <- spm_continuous(dat=data,a=-0.05, f1=80, Q=2e-8, f=80, b=5, mu0=2e-5, theta=0.08)
## Parameter a achieved lower/upper bound.
## -0.055
## Parameter Q achieved lower/upper bound.
## 2.2e-08
## Parameter f achieved lower/upper bound.
## 88
## Parameter mu0 achieved lower/upper bound.
## 2.2e-05
pars
## $a ## [,1] ## [1,] -0.055 ## ##$f1
##             [,1]
## [1,] 79.61017189
##
## $Q ## [,1] ## [1,] 2.2e-08 ## ##$f
##      [,1]
## [1,]   88
##
## $b ## [,1] ## [1,] 5.03173818 ## ##$mu0
## [1] 2.2e-05
##
## $theta ## [1] 0.07951200763 ## ##$status
## [1] 5
##
## $LogLik ## [1] -7080.5778 ## ##$objective
## [1] 7080.572992
##
## $message ## [1] "NLOPT_MAXEVAL_REACHED: Optimization stopped because maxeval (above) was reached." ## ##$limit
## [1] TRUE
##
## attr(,"class")
## [1] "spm.continuous"

### Coefficient conversion between continuous- and discrete-time models

The coefficient conversion between continuous- and discrete-time models is as follows (‘c’ and ‘d’ denote continuous- and discrete-time models respectively; note: these equations can be used if intervals between consecutive observations of discrete- and continuous-time models are equal; it also required that matrices $$\mathbf{a}_c$$ and $$\mathbf{Q}_{c,d}$$ must be full-rank matrices):

$$\mathbf{Q}_c = \mathbf{Q}_d$$

$$\mathbf{a}_c = \mathbf{R}_d - I(k)$$

$$\mathbf{b}_c = \mathbf{\Sigma}$$

$${\mathbf{f}_1}_c = -\mathbf{a}_c^{-1} \times \mathbf{u}_d$$

$$\mathbf{f}_c = -0.5 \mathbf{b}_d \times \mathbf{Q}^{-1}_d$$

$${\mu_0}_c = {\mu _0}_d - \mathbf{f}_c \times \mathbf{Q_c} \times \mathbf{f}_c^*$$

$$\theta_c = \theta_d$$

where $$k$$ is a number of covariates, which is equal to model’s dimension and ’*’ denotes transpose operation; $$\mathbf{\Sigma}$$ is a $$k \times 1$$ matrix which contains s.d.s of corresponding residuals (residuals of a linear regression $$\mathbf{Y}(t+1) = \mathbf{u} + \mathbf{R}\mathbf{Y}(t) + \mathbf{\epsilon}$$; s.d. is a standard deviation), $$I(k)$$ is an identity $$k \times k$$ matrix.

### Model with time-dependent coefficients

In previous models, we assumed that coefficients is sort of time-dependant: we multiplied them on to $$e^{\theta t}$$. In general, this may not be the case (A. I. Yashin, Arbeev, Kulminski, et al. 2007). We extend this to a general case, i.e. (we consider one-dimensional case):

$$\mathbf{a(t)} = \mathbf{par}_1 t + \mathbf{par}_2$$ - linear function.

The corresponding equations will be equivalent to one-dimensional continuous case described above.

#### Example

library(stpm)
#Data preparation:
n <- 10
data <- simdata_time_dep(N=n)
# Estimation:
opt.par <- spm_time_dep(data,
start = list(a = -0.05, f1 = 80, Q = 2e-08, f = 80, b = 5, mu0 = 0.001),
frm = list(at = "a", f1t = "f1", Qt = "Q", ft = "f", bt = "b", mu0t= "mu0"))
opt.par
## [[1]]
## [[1]]$a ## [1] -0.05 ## ## [[1]]$f1
## [1] 80
##
## [[1]]$Q ## [1] 2e-08 ## ## [[1]]$f
## [1] 80
##
## [[1]]$b ## [1] 5 ## ## [[1]]$mu0
## [1] 0.001

## Setting lower and upper boundaries of the model parameters

Lower and upper boundaries can be set up with parameters $$lb$$ and $$ub$$, which represents simple numeric vectors. Note: lengths of $$lb$$ and $$ub$$ must be the same as the total length of the parameters. Lower and upper boundaries can be set for continuous-time and time-dependent models only.

### Setting lb and ub for continuous-time model

#### One covariate

Below we show the example of setting up $$lb$$ and $$ub$$ when we have a single covariate:

library(stpm)
data <- simdata_cont(N=10, ystart = 80, a = -0.1, Q = 1e-06, mu0 = 1e-5, theta = 0.08, f1 = 80, f=80, b=1, dt=1, sd0=5)
ans <- spm_continuous(dat=data,
a = -0.1,
f1 = 82,
Q = 1.4e-6,
f = 77,
b = 1,
mu0 = 1.6e-5,
theta = 0.1,
stopifbound = FALSE,
lb=c(-0.2, 60, 0.1e-6, 60, 0.1, 0.1e-5, 0.01),
ub=c(0, 140, 5e-06, 140, 3, 5e-5, 0.20))
ans
## $a ## [,1] ## [1,] -0.09498756473 ## ##$f1
##             [,1]
## [1,] 80.59501818
##
## $Q ## [,1] ## [1,] 3.512785119e-06 ## ##$f
##             [,1]
## [1,] 117.0678995
##
## $b ## [,1] ## [1,] 1.018411166 ## ##$mu0
## [1] 4.254766784e-05
##
## $theta ## [1] 0.1887345857 ## ##$status
## [1] 5
##
## $LogLik ## [1] -685.7213874 ## ##$objective
## [1] 685.6869249
##
## $message ## [1] "NLOPT_MAXEVAL_REACHED: Optimization stopped because maxeval (above) was reached." ## ##$limit
## [1] FALSE
##
## attr(,"class")
## [1] "spm.continuous"

#### Two covariates

This is an example for two physiological variables (covariates).

library(stpm)

data <- simdata_cont(N=10,
a=matrix(c(-0.1,  0.001, 0.001, -0.1), nrow = 2, ncol = 2, byrow = T),
f1=t(matrix(c(100, 200), nrow = 2, ncol = 1, byrow = F)),
Q=matrix(c(1e-06, 1e-7, 1e-7,  1e-06), nrow = 2, ncol = 2, byrow = T),
f=t(matrix(c(100, 200), nrow = 2, ncol = 1, byrow = F)),
b=matrix(c(1, 2), nrow = 2, ncol = 1, byrow = F),
mu0=1e-4,
theta=0.08,
ystart = c(100,200), sd0=c(5, 10), dt=1)

a.d <- matrix(c(-0.15,  0.002, 0.002, -0.15), nrow = 2, ncol = 2, byrow = T)
f1.d <- t(matrix(c(95, 195), nrow = 2, ncol = 1, byrow = F))
Q.d <- matrix(c(1.2e-06, 1.2e-7, 1.2e-7,  1.2e-06), nrow = 2, ncol = 2, byrow = T)
f.d <- t(matrix(c(105, 205), nrow = 2, ncol = 1, byrow = F))
b.d <- matrix(c(1, 2), nrow = 2, ncol = 1, byrow = F)
mu0.d <- 1.1e-4
theta.d <- 0.07

ans <- spm_continuous(dat=data,
a = a.d,
f1 = f1.d,
Q = Q.d,
f = f.d,
b = b.d,
mu0 = mu0.d,
theta = theta.d,
lb=c(-0.5, ifelse(a.d[2,1] > 0, a.d[2,1]-0.5*a.d[2,1], a.d[2,1]+0.5*a.d[2,1]), ifelse(a.d[1,2] > 0, a.d[1,2]-0.5*a.d[1,2], a.d[1,2]+0.5*a.d[1,2]), -0.5,
80, 100,
Q.d[1,1]-0.5*Q.d[1,1], ifelse(Q.d[2,1] > 0, Q.d[2,1]-0.5*Q.d[2,1], Q.d[2,1]+0.5*Q.d[2,1]), ifelse(Q.d[1,2] > 0, Q.d[1,2]-0.5*Q.d[1,2], Q.d[1,2]+0.5*Q.d[1,2]), Q.d[2,2]-0.5*Q.d[2,2],
80, 100,
0.1, 0.5,
0.1e-4,
0.01),
ub=c(-0.08,  0.002,  0.002, -0.08,
110, 220,
Q.d[1,1]+0.1*Q.d[1,1], ifelse(Q.d[2,1] > 0, Q.d[2,1]+0.1*Q.d[2,1], Q.d[2,1]-0.1*Q.d[2,1]), ifelse(Q.d[1,2] > 0, Q.d[1,2]+0.1*Q.d[1,2], Q.d[1,2]-0.1*Q.d[1,2]), Q.d[2,2]+0.1*Q.d[2,2],
110, 220,
1.5, 2.5,
1.2e-4,
0.10))
ans
## $a ## [,1] [,2] ## [1,] -0.150605482175 0.001679541154 ## [2,] 0.001975476789 -0.148864241454 ## ##$f1
##             [,1]
## [1,] 106.0260503
## [2,] 193.4146154
##
## $Q ## [,1] [,2] ## [1,] 1.307268804e-06 1.303231265e-07 ## [2,] 1.320000000e-07 1.320000000e-06 ## ##$f
##             [,1]
## [1,] 109.8219998
## [2,] 216.8226510
##
## $b ## [,1] ## [1,] 1.251339994 ## [2,] 1.933924800 ## ##$mu0
## [1] 0.0001141832977
##
## $theta ## [1] 0.07653264804 ## ##$status
## [1] 5
##
## $LogLik ## [1] 1565.568825 ## ##$objective
## [1] -2229.740166
##
## $message ## [1] "NLOPT_MAXEVAL_REACHED: Optimization stopped because maxeval (above) was reached." ## ##$limit
## [1] FALSE
##
## attr(,"class")
## [1] "spm.continuous"

### Setting lb and ub for model with time-dependent coefficients

This model uses only one covariate, therefore setting-up model parameters is easy:

n <- 10
data <- simdata_time_dep(N=n)
# Estimation:
opt.par <- spm_time_dep(data, start=list(a=-0.05, f1=80, Q=2e-08, f=80, b=5, mu0=0.001),
lb=c(-1, 30, 1e-8, 30, 1, 1e-6), ub=c(0, 120, 5e-8, 130, 10, 1e-2))
opt.par
## [[1]]
## [[1]]$a ## [1] -0.05 ## ## [[1]]$f1
## [1] 80
##
## [[1]]$Q ## [1] 2e-08 ## ## [[1]]$f
## [1] 80
##
## [[1]]$b ## [1] 5 ## ## [[1]]$mu0
## [1] 0.001

#### Special case when some model parameter functions are equal to zero

Imagine a situation when one parameter function you want to be equal to zero: $$f=0$$. Let’s emulate this case:

library(stpm)
n <- 10
data <- simdata_time_dep(N=n)
# Estimation:
opt.par <- spm_time_dep(data, frm = list(at="a", f1t="f1", Qt="Q", ft="0", bt="b", mu0t="mu0"))
opt.par
## [[1]]
## [[1]]$a ## [1] -0.05 ## ## [[1]]$f1
## [1] 80
##
## [[1]]$Q ## [1] 2e-08 ## ## [[1]]$b
## [1] 80
##
## [[1]]$mu0 ## [1] 5 As you can see, there is no parameter $$f$$ in $$opt.par$$. This because we set $$f=0$$ in $$frm$$! Then, is you want to set the constraints, you must not specify the starting value (parameter $$start$$) and $$lb$$/$$ub$$ for the parameter $$f$$ (otherwise, the function raises an error): n <- 10 data <- simdata_time_dep(N=n) # Estimation: opt.par <- spm_time_dep(data, frm = list(at="a", f1t="f1", Qt="Q", ft="0", bt="b", mu0t="mu0"), start=list(a=-0.05, f1=80, Q=2e-08, b=5, mu0=0.001), lb=c(-1, 30, 1e-8, 1, 1e-6), ub=c(0, 120, 5e-8, 10, 1e-2)) opt.par ## [[1]] ## [[1]]$a
## [1] -0.05
##
## [[1]]$f1 ## [1] 80 ## ## [[1]]$Q
## [1] 2e-08
##
## [[1]]$b ## [1] 5 ## ## [[1]]$mu0
## [1] 0.001

You can do the same manner if you want two or more parameters to be equal to zero.

## Simulation (individual trajectory projection, also known as microsimulations)

We added one- and multi- dimensional simulation to be able to generate test data for hyphotesis testing. Data, which can be simulated can be discrete (equal intervals between observations) and continuous (with arbitrary intervals).

### Discrete-time simulation

The corresponding function is (k - a number of variables(covariates), equal to model’s dimension):

simdata_discr(N=100, a=-0.05, f1=80, Q=2e-8, f=80, b=5, mu0=1e-5, theta=0.08, ystart=80, tstart=30, tend=105, dt=1)

Here:

N - Number of individuals

a - A matrix of kxk, which characterize the rate of the adaptive response

f1 - A particular state, which if a deviation from the normal (or optimal). This is a vector with length of k

Q - A matrix of k by k, which is a non-negative-definite symmetric matrix

f - A vector-function (with length k) of the normal (or optimal) state

b - A diffusion coefficient, k by k matrix

mu0 - mortality at start period of time (baseline hazard)

theta - A displacement coefficient of the Gompertz function

ystart - A vector with length equal to number of dimensions used, defines starting values of covariates

tstart - A number that defines a start time (30 by default). Can be a number (30 by default) or a vector of two numbers: c(a, b) - in this case, starting value of time is simulated via uniform(a,b) distribution.

tend - A number, defines a final time (105 by default)

dt - A time interval between observations.

This function returns a table with simulated data, as shown in example below:

library(stpm)
data <- simdata_discr(N=10)
head(data)
##   id xi t1 t2          y1     y1.next
## 1  1  0 30 31 77.96865128 74.95020713
## 2  1  0 31 32 74.95020713 78.81422683
## 3  1  0 32 33 78.81422683 77.18093430
## 4  1  0 33 34 77.18093430 78.93473268
## 5  1  0 34 35 78.93473268 82.75382076
## 6  1  0 35 36 82.75382076 77.95326935

### Continuous-time simulation

The corresponding function is (k - a number of variables(covariates), equal to model’s dimension):

simdata_cont(N=100, a=-0.05, f1=80, Q=2e-07, f=80, b=5, mu0=2e-05, theta=0.08, ystart=80, tstart=c(30,50), tend=105)

Here:

N - Number of individuals

a - A matrix of kxk, which characterize the rate of the adaptive response

f1 - A particular state, which if a deviation from the normal (or optimal). This is a vector with length of k

Q - A matrix of k by k, which is a non-negative-definite symmetric matrix

f - A vector-function (with length k) of the normal (or optimal) state

b - A diffusion coefficient, k by k matrix

mu0 - mortality at start period of time (baseline hazard)

theta - A displacement coefficient of the Gompertz function

ystart - A vector with length equal to number of dimensions used, defines starting values of covariates

tstart - A number that defines a start time (30 by default). Can be a number (30 by default) or a vector of two numbers: c(a, b) - in this case, starting value of time is simulated via uniform(a,b) distribution.

tend - A number, defines a final time (105 by default)

This function returns a table with simulated data, as shown in example below:

library(stpm)
data <- simdata_cont(N=10)
head(data)
##   id xi          t1          t2          y1     y1.next
## 1  0  0 32.91452991 34.59870981 80.33014489 79.31869904
## 2  0  0 34.59870981 36.42590572 79.31869904 86.42858780
## 3  0  0 36.42590572 37.64290358 86.42858780 86.12034765
## 4  0  0 37.64290358 38.85700475 86.12034765 82.89149789
## 5  0  0 38.85700475 40.69648447 82.89149789 74.78112714
## 6  0  0 40.69648447 42.64176469 74.78112714 69.82268255

## SPM with partially observed covariates

Stochastic Process Model has many applications in analysis of longitudinal biodemographic data. Such data contain various physiological variables (known as covariates). Data can also potentially contain genetic information available for all or a part of participants. Taking advantage from both genetic and non-genetic information can provide future insights into a broad range of processes describing aging-related changes in the organism.

### Method

In this package, SPM with partially observed covariates is implemented in form of GenSPM (Genetic SPM), presented in (Arbeev et al. 2009) and further advanced in (Arbeev et al. 2014), further elaborates the basic stochastic process model conception by introducing a categorical variable, $$Z$$, which may be a specific value of a genetic marker or, in general, any categorical variable. Currently, $$Z$$ has two gradations: 0 or 1 in a genetic group of interest, assuming that $$P(Z=1) = p$$, $$p \in [0, 1]$$, were $$p$$ is the proportion of carriers and non-carriers of an allele in a population. Example of longitudinal data with genetic component $$Z$$ is provided below.

library(stpm)
data <- sim_pobs(N=10)
head(data)
##   id xi          t1          t2 Z          y1     y1.next
## 1  0  0 81.70579868 82.63613491 0 78.34545722 77.70714589
## 2  0  0 82.63613491 83.64611170 0 77.70714589 81.26466661
## 3  0  0 83.64611170 84.58459104 0 81.26466661 88.79834425
## 4  0  0 84.58459104 85.56018957 0 88.79834425 90.16389169
## 5  0  0 85.56018957 86.46097239 0 90.16389169 91.15936311
## 6  0  0 86.46097239 87.48069276 0 91.15936311 86.35212341

In the specification of the SPM described in 2007 paper by Yashin and colleagues (A. I. Yashin, Arbeev, Akushevich, et al. 2007) the stochastic differential equation describing the age dynamics of a physiological variable (a dynamic component of the model) is:

$$dY(t) = a(Z, t)(Y(t) - f1(Z, t))dt + b(Z, t)dW(t), Y(t = t_0)$$

Here in this equation, $$Y(t)$$ is a $$k \times 1$$ matrix, where $$k$$ is a number of covariates, which is a model dimension) describing the value of a physiological variable at a time (e.g. age) t. $$f_1(Z,t)$$ is a $$k \times 1$$ matrix that corresponds to the long-term average value of the stochastic process $$Y(t)$$, which describes a trajectory of individual variable influenced by different factors represented by a random Wiener process $$W(t)$$. The negative feedback coefficient $$a(Z,t)$$ ($$k \times k$$ matrix) characterizes the rate at which the stochastic process goes to its mean. In research on aging and well-being, $$f_1(Z,t)$$ represents the average allostatic trajectory and $$a(t)$$ in this case represents the adaptive capacity of the organism. Coefficient $$b(Z,t)$$ ($$k \times 1$$ matrix) characterizes a strength of the random disturbances from Wiener process $$W(t)$$. All of these parameters depend on $$Z$$ (a genetic marker having values 1 or 0). The following function $$\mu(t,Y(t))$$ represents a hazard rate:

$$\mu(t,Y(t)) = \mu_0(t) + (Y(t) - f(Z, t))^*Q(Z, t)(Y(t) - f(Z, t))$$

In this equation: $$\mu_0(t)$$ is the baseline hazard, which represents a risk when $$Y(t)$$ follows its optimal trajectory; f(t) ($$k \times 1$$ matrix) represents the optimal trajectory that minimizes the risk and $$Q(Z, t)$$ ($$k \times k$$ matrix) represents a sensitivity of risk function to deviation from the norm. In general, model coefficients $$a(Z, t)$$, $$f1(Z, t)$$, $$Q(Z, t)$$, $$f(Z, t)$$, $$b(Z, t)$$ and $$\mu_0(t)$$ are time(age)-dependent. Once we have data, we then can run analysis, i.e. estimate coefficients (they are assumed to be time-independent and data here is simulated):

library(stpm)
#Generating data:
data <- sim_pobs(N=10)
head(data)
##   id xi          t1          t2 Z          y1     y1.next
## 1  0  0 73.12292986 74.06692660 0 80.70401805 88.25724510
## 2  0  0 74.06692660 75.01785380 0 88.25724510 84.09597770
## 3  0  0 75.01785380 76.06218815 0 84.09597770 78.73632764
## 4  0  0 76.06218815 77.05379218 0 78.73632764 79.45687786
## 5  0  0 77.05379218 78.02605115 0 79.45687786 78.93149440
## 6  0  0 78.02605115 79.01875033 0 78.93149440 69.89422735
#Parameters estimation:
pars <- spm_pobs(x=data)
pars
## $aH ## [,1] ## [1,] -0.05130364936 ## ##$aL
##                [,1]
## [1,] -0.01024774755
##
## $f1H ## [,1] ## [1,] 61.5318577 ## ##$f1L
##             [,1]
## [1,] 72.04588562
##
## $QH ## [,1] ## [1,] 2.109557233e-08 ## ##$QL
##                 [,1]
## [1,] 2.262883621e-08
##
## $fH ## [,1] ## [1,] 65.89856393 ## ##$fL
##             [,1]
## [1,] 87.29866489
##
## $bH ## [,1] ## [1,] 3.632520911 ## ##$bL
##             [,1]
## [1,] 5.046458225
##
## $mu0H ## [1] 8.790171786e-06 ## ##$mu0L
## [1] 9.003907871e-06
##
## $thetaH ## [1] 0.072213777 ## ##$thetaL
## [1] 0.09000393391
##
## $p ## [1] 0.255052107 ## ##$limit
## [1] FALSE
##
## attr(,"class")
## [1] "pobs.spm"

Here and represents parameters when $$Z$$ = 1 (H) and 0 (L).

### Joint analysis of two datasets: first dataset with genetic and second dataset with non-genetic component

library(stpm)
data.genetic <- sim_pobs(N=5, mode='observed')
head(data.genetic)
##   id xi          t1          t2 Z          y1     y1.next
## 1  0  0 81.30261515 82.37815787 0 81.86051501 77.36316584
## 2  0  0 82.37815787 83.37960803 0 77.36316584 72.99358524
## 3  0  0 83.37960803 84.31653383 0 72.99358524 68.47142336
## 4  0  0 84.31653383 85.22570729 0 68.47142336 69.05564899
## 5  0  0 85.22570729 86.20435298 0 69.05564899 64.38584330
## 6  0  0 86.20435298 87.20720174 0 64.38584330 64.17757900
data.nongenetic <- sim_pobs(N=10, mode='unobserved')
head(data.nongenetic)
##   id xi           t1           t2          y1     y1.next
## 1  0  0 101.40351475 102.43415759 78.73900402 70.84733959
## 2  0  0 102.43415759 103.38566410 70.84733959 68.82386758
## 3  0  0 103.38566410 104.36668151 68.82386758 69.07764575
## 4  1  0  85.42064926  86.49500468 79.76337453 82.22952762
## 5  1  0  86.49500468  87.49709424 82.22952762 72.33267890
## 6  1  0  87.49709424  88.50258890 72.33267890 75.98092900
#Parameters estimation:
pars <- spm_pobs(x=data.genetic, y = data.nongenetic, mode='combined')
pars
## $aH ## [,1] ## [1,] -0.05436391599 ## ##$aL
##                [,1]
## [1,] -0.01082216864
##
## $f1H ## [,1] ## [1,] 54.14953148 ## ##$f1L
##             [,1]
## [1,] 79.58470122
##
## $QH ## [,1] ## [1,] 1.842585008e-08 ## ##$QL
##                 [,1]
## [1,] 1.998741701e-08
##
## $fH ## [,1] ## [1,] 65.46246282 ## ##$fL
##            [,1]
## [1,] 86.3317748
##
## $bH ## [,1] ## [1,] 3.924976631 ## ##$bL
##           [,1]
## [1,] 4.9760502
##
## $mu0H ## [1] 8.780659142e-06 ## ##$mu0L
## [1] 9.03116267e-06
##
## $thetaH ## [1] 0.07228407181 ## ##$thetaL
## [1] 0.09016449506
##
## $p ## [1] 0.2729964624 ## ##$limit
## [1] FALSE
##
## attr(,"class")
## [1] "pobs.spm"

Here mode ‘observed’ is used for simlation of data with genetic component $$Z$$ and ‘unobserved’ - without genetic component.

## Genetic SPM ‘GSPM’

This type of SPM also uses genetic component by analogy from the previous chapters but uses explicit gradient function which speeds up computations significantly. See (He et al. 2017) for details. Below we provide examples of usage:

library(stpm)
data(ex_spmcon1dg)
head(ex_data$spm_data) ## id xi t1 t2 y y.next ## 1 1 0 30 31 2.000000000 2.024328135 ## 2 1 0 31 32 2.024328135 1.927486318 ## 3 1 0 32 33 1.927486318 1.899083801 ## 4 1 0 33 34 1.899083801 2.061574385 ## 5 1 0 34 35 2.061574385 2.034558435 ## 6 1 0 35 36 2.034558435 2.114382051 head(ex_data$gene_data)
##   id geno
## 1  1    1
## 2  2    1
## 3  3    0
## 4  4    0
## 5  5    1
## 6  6    0
res <- spm_con_1d_g(spm_data=ex_data$spm_data, gene_data=ex_data$gene_data,
a = -0.02, b=0.2, q=0.01, f=3, f1=3, mu0=0.01, theta=1e-05,
upper=c(-0.01,3,0.1,10,10,0.1,1e-05), lower=c(-1,0.01,0.00001,1,1,0.001,1e-05),
effect=c('q'))
## [1] "Initial values:"
##  [1] -2e-02 -2e-02  2e-01  2e-01  1e-02  1e-02  3e+00  3e+00  3e+00  1e-02
## [11]  1e-02  1e-05
## [1] "Lower bounds:"
##  [1] -1e+00 -1e+00  1e-02  1e-02  1e-05  1e-05  1e+00  1e+00  1e+00  1e-03
## [11]  1e-03  1e-05
## [1] "Upper bounds:"
##  [1] -1e-02 -1e-02  3e+00  3e+00  1e-01  1e-01  1e+01  1e+01  1e+01  1e-01
## [11]  1e-01  1e-05
res
## $est ## a b q_0 q_2 f ## [1,] -0.03054228766 0.101197293 0.00448415079 0.004879437494 2.045171446 ## f1 mu0 theta ## [1,] 3.010584769 0.001420354314 1e-05 ## ##$lik
## [1] -121717.8478
##
## $con ## [1] 3 ## ##$message
## [1] "NLOPT_FTOL_REACHED: Optimization stopped because ftol_rel or ftol_abs (above) was reached."
##
## $hessian ## [,1] [,2] [,3] [,4] ## [1,] 2418684.83379972 6.662849963e+05 1986.7582655 1121.9246461 ## [2,] 666284.99625638 1.332736169e+07 -723.6350768 -551.1722288 ## [3,] 1986.75826554 -7.236350768e+02 2667427.4599746 855374.8268798 ## [4,] 1121.92464606 -5.511722288e+02 855374.8268798 855374.8268798 ## [5,] -82.61078654 -3.786615040e+01 -32507.4291866 -16187.8855511 ## [6,] -85019.08324387 -9.944962282e+01 -247.8411941 -131.6918542 ## [7,] -11704.11854947 1.410671930e+04 4620297.2574421 2397280.1060920 ## [8,] -895.86209312 1.023340269e+03 468707.7861779 232112.4802746 ## [,5] [,6] [,7] [,8] ## [1,] -82.610786536 -85019.083243867 -11704.118549 -895.8620931 ## [2,] -37.866150397 -99.449622818 14106.719296 1023.3402687 ## [3,] -32507.429186647 -247.841194103 4620297.257442 468707.7861779 ## [4,] -16187.885551148 -131.691854199 2397280.106092 232112.4802746 ## [5,] 551.678053337 1.016137425 -90004.673127 -8475.9149129 ## [6,] 1.016137425 6201.557338306 1347.685889 126.8541484 ## [7,] -90004.673126606 1347.685889021 38250575.535221 2134987.4964705 ## [8,] -8475.914912931 126.854148429 2134987.496470 215962.9544164 ## ##$beta
##      beta_a beta_b          beta_q beta_f beta_mu0
## [1,]     NA     NA 0.0001976433519     NA       NA

Here: spm_data - A dataset for the SPM model. See the STPM package for more details about the format.

gene_data - A two column dataset containing the genotypes for the individuals in spm_data. The first column id is the ID of the individuals in dataset spm_data, and the second column geno is the genotype.

a - The initial value for the paramter . The initial value will be predicted if not specified.

b - The initial value for the paramter . The initial value will be predicted if not specified.

q - The initial value for the paramter . The initial value will be predicted if not specified.

f - The initial value for the paramter . The initial value will be predicted if not specified.

f1 - The initial value for the paramter . The initial value will be predicted if not specified.

mu0 - The initial value for the paramter in the baseline hazard. The initial value will be predicted if not specified.

theta - The initial value for the paramter in the baseline hazard. The initial value will be predicted if not specified.

lower - A vector of the lower bound of the parameters.

upper - A vector of the upper bound of the parameters.

effect - A character vector of the parameters that are linked to genotypes. The vector can contain any combination of , , , , .

control - A list of the control parameters for the optimization paramters.

global - A logical variable indicating whether the MLSL (TRUE) or the L-BFGS (FALSE) algorithm is used for the optimization.

verbose - A logical variable indicating whether initial information is printed.

ahessian - A logical variable indicating whether the approximate (FALSE) or analytical (TRUE) Hessian is returned.

est - The estimates of the parameters.

hessian - The Hessian matrix of the estimates.

lik - The minus log-likelihood.

con - A number indicating the convergence. See the ‘nloptr’ package for more details.

message - Extra message about the convergence. See the ‘nloptr’ package for more details.

beta - The coefficients of the genetic effect on the parameters to be linked to genotypes.

## Multiple imputation with spm.impute(…)

The SPM offers longitudinal data imputation with results that are better than from other imputation tools since it preserves data structure, i.e. relation between Y(t) and mu(Y(t),t). Below there are two examples of multiple data imputation with function spm.impute(…).

library(stpm)

#######################################################################
############## One dimensional case (one covariate) ###################
#######################################################################

## Data preparation (short format)#
data <- simdata_discr(N=1000, dt = 2, format="short")

miss.id <- sample(x=dim(data)[1], size=round(dim(data)[1]/4)) # ~25% missing data
incomplete.data <- data
incomplete.data[miss.id,4] <- NA
# End of data preparation #

##### Multiple imputation with SPM #####
imp.data <- spm.impute(x=incomplete.data, id=1, case="xi", t1=3, covariates="y1", minp=1, theta_range=seq(0.075, 0.09, by=0.001))$imputed ##### Look at the incomplete data with missings ##### head(incomplete.data) ## id xi t y1 ## 1 1 0 30 79.30322156 ## 2 1 0 32 NA ## 3 1 0 34 80.97494790 ## 4 1 0 36 NA ## 5 1 0 38 65.22183053 ## 6 1 0 40 74.79238922 ##### Look at the imputed data ##### head(imp.data) ## id xi t y1 ## 1 1 0 30 79.30322156 ## 2 1 0 32 79.32755390 ## 3 1 0 34 80.97494790 ## 4 1 0 36 80.92187520 ## 5 1 0 38 65.22183053 ## 6 1 0 40 74.79238922 ######################################################### ################ Two-dimensional case ################### ######################################################### ## Parameters for data simulation # a <- matrix(c(-0.05, 0.01, 0.01, -0.05), nrow=2) f1 <- matrix(c(90, 30), nrow=1, byrow=FALSE) Q <- matrix(c(1e-7, 1e-8, 1e-8, 1e-7), nrow=2) f0 <- matrix(c(80, 25), nrow=1, byrow=FALSE) b <- matrix(c(5, 3), nrow=2, byrow=TRUE) mu0 <- 1e-04 theta <- 0.07 ystart <- matrix(c(80, 25), nrow=2, byrow=TRUE) ## Data preparation # data <- simdata_discr(N=1000, a=a, f1=f1, Q=Q, f=f0, b=b, ystart=ystart, mu0 = mu0, theta=theta, dt=2, format="short") ## Delete some observations in order to have approx. 25% missing data incomplete.data <- data miss.id <- sample(x=dim(data)[1], size=round(dim(data)[1]/4)) incomplete.data <- data incomplete.data[miss.id,4] <- NA miss.id <- sample(x=dim(data)[1], size=round(dim(data)[1]/4)) incomplete.data[miss.id,5] <- NA ## End of data preparation # ###### Multiple imputation with SPM ##### imp.data <- spm.impute(x=incomplete.data, id=1, case="xi", t1=3, covariates=c("y1", "y2"), minp=1, theta_range=seq(0.060, 0.07, by=0.001))$imputed

###### Look at the incomplete data with missings #####
head(incomplete.data)
##   id xi  t          y1          y2
## 1  1  0 30 79.88608631 22.01654900
## 2  1  0 32          NA 21.91146314
## 3  1  0 34 68.35619345          NA
## 4  1  0 36 65.85888119 21.96011119
## 5  1  0 38 67.34629181 18.93992585
## 6  1  0 40          NA 24.45030784
###### Look at the imputed data #####
head(imp.data)
##   id xi  t          y1          y2
## 1  1  0 30 79.88608631 22.01654900
## 2  1  0 32 80.18366311 21.91146314
## 3  1  0 34 68.35619345 22.78051949
## 4  1  0 36 65.85888119 21.96011119
## 5  1  0 38 67.34629181 18.93992585
## 6  1  0 40 68.22541009 24.45030784

## References

Akushevich, I., A. Kulminski, and K. G. Manton. 2005. “Life Tables with Covariates: Dynamic Model for Nonlinear Analysis of Longitudinal Data.” Mathematical Population Studies 12 (2). Informa UK Limited: 51–80. doi:10.1080/08898480590932296.

Arbeev, Konstantin G., Igor Akushevich, Alexander M. Kulminski, Svetlana V. Ukraintseva, and Anatoliy I. Yashin. 2014. “Joint Analyses of Longitudinal and Time-to-Event Data in Research on Aging: Implications for Predicting Health and Survival.” Frontiers in Public Health 2 (November). Frontiers Media SA. doi:10.3389/fpubh.2014.00228.

Arbeev, Konstantin G., Igor Akushevich, Alexander M. Kulminski, Liubov S. Arbeeva, Lucy Akushevich, Svetlana V. Ukraintseva, Irina V. Culminskaya, and Anatoli I. Yashin. 2009. “Genetic Model for Longitudinal Studies of Aging, Health, and Longevity and Its Potential Application to Incomplete Data.” Journal of Theoretical Biology 258 (1). Elsevier BV: 103–11. doi:10.1016/j.jtbi.2009.01.023.

He, Liang, Ilya Zhbannikov, Konstantin G. Arbeev, Anatoliy I. Yashin, and Alexander M. Kulminski. 2017. “A Genetic Stochastic Process Model for Genome-Wide Joint Analysis of Biomarker Dynamics and Disease Susceptibility with Longitudinal Data.” Genetic Epidemiology 41 (7). Wiley: 620–35. doi:10.1002/gepi.22058.

Woodbury, Max A., and Kenneth G. Manton. 1977. “A Random-Walk Model of Human Mortality and Aging.” Theoretical Population Biology 11 (1). Elsevier BV: 37–48. doi:10.1016/0040-5809(77)90005-3.

Yashin, Anatoli I., Konstantin G. Arbeev, Igor Akushevich, Aliaksandr Kulminski, Lucy Akushevich, and Svetlana V. Ukraintseva. 2007. “Stochastic Model for Analysis of Longitudinal Data on Aging and Mortality.” Mathematical Biosciences 208 (2). Elsevier BV: 538–51. doi:10.1016/j.mbs.2006.11.006.

Yashin, Anatoli I., Konstantin G. Arbeev, Aliaksandr Kulminski, Igor Akushevich, Lucy Akushevich, and Svetlana V. Ukraintseva. 2007. “Health Decline, Aging and Mortality: How Are They Related?” Biogerontology 8 (3). Springer Nature: 291–302. doi:10.1007/s10522-006-9073-3.