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.
install.packages("stpm")require(devtools)
devtools::install_github("izhbannikov/stpm")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.89139The packate accepts longitudinal data in two formats: “short” and “long”.
## id xi t y1
## 1 1 0 30 81.23557
## 2 1 0 31 82.28856
## 3 1 0 32 87.25457
## 4 1 0 33 87.40676
## 5 1 0 34 87.12992
## 6 1 0 35 86.49620## id xi t1 t2 y1 y1.next
## 1 1 0 30 31 84.11381 84.54564
## 2 1 0 31 32 84.54564 93.94254
## 3 1 0 32 33 93.94254 94.66992
## 4 1 0 33 34 94.66992 89.04257
## 5 1 0 34 35 89.04257 90.19970
## 6 1 0 35 36 90.19970 86.31144There 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.78035 96.13678
## 2 1 0 31 32 96.13678 93.06003
## 3 1 0 32 33 93.06003 93.77237
## 4 1 0 33 34 93.77237 102.77195
## 5 1 0 34 35 102.77195 105.63940
## 6 1 0 35 36 105.63940 100.54099In 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 36.45512 80.21420
## 2 0 0 38.28591 85.94959
## 3 0 0 39.88678 86.22148
## 4 0 0 41.42762 84.12783
## 5 0 0 42.73823 84.71374
## 6 0 0 44.58876 93.88350The 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}\)).
library(stpm)
#Data simulation (200 individuals)
data <- simdata_discr(N=100)
#Estimation of parameters
pars <- spm_discrete(data)
pars## $dmodel
## $dmodel$theta
## [1] 0.086
##
## $dmodel$mu0
## [1] 6.624258952e-05
##
## $dmodel$b
## [1] -1.411133262e-06
##
## $dmodel$Q
## [,1]
## [1,] 8.553331358e-09
##
## $dmodel$u
## [1] 3.66334276
##
## $dmodel$u.std.err
## (Intercept)
## 0.3256923735
##
## $dmodel$R
## [,1]
## [1,] 0.9541720023
##
## $dmodel$R.std.err
## y1_1
## [1,] 0.00400638158
##
## $dmodel$Sigma
## [1] 4.996513374
##
##
## $cmodel
## $cmodel$a
## [,1]
## [1,] -0.04582799769
##
## $cmodel$f1
## [,1]
## [1,] 79.93678417
##
## $cmodel$Q
## [,1]
## [1,] 8.553331358e-09
##
## $cmodel$f
## [,1]
## [1,] 82.49027212
##
## $cmodel$b
## [,1]
## [1,] 4.996513374
##
## $cmodel$mu0
## [,1]
## [1,] 8.040206117e-06
##
## $cmodel$theta
## [1] 0.086
##
##
## attr(,"class")
## [1] "spm.discrete"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.
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 31.54041716 32.94180463 81.16635305 78.63119219
## 2 0 0 32.94180463 34.72039120 78.63119219 79.48928398
## 3 0 0 34.72039120 36.68805054 79.48928398 79.98422050
## 4 0 0 36.68805054 38.27744815 79.98422050 84.10036985
## 5 0 0 38.27744815 40.13373057 84.10036985 81.88510408
## 6 0 0 40.13373057 41.51573423 81.88510408 90.17307264#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)
pars## $a
## [,1]
## [1,] -0.05499599743
##
## $f1
## [,1]
## [1,] 80.85351159
##
## $Q
## [,1]
## [1,] 2.199853058e-08
##
## $f
## [,1]
## [1,] 87.99486979
##
## $b
## [,1]
## [1,] 4.995789558
##
## $mu0
## [1] 2.2e-05
##
## $theta
## [1] 0.08252151133
##
## $status
## [1] 5
##
## $LogLik
## [1] -6589.039843
##
## $objective
## [1] 6589.039634
##
## $message
## [1] "NLOPT_MAXEVAL_REACHED: Optimization stopped because maxeval (above) was reached."
##
## $limit
## [1] FALSE
##
## attr(,"class")
## [1] "spm.continuous"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.
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.
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## $a
## [1] -0.05
##
## $f1
## [1] 80
##
## $Q
## [1] 2e-08
##
## $f
## [1] 80
##
## $b
## [1] 5
##
## $mu0
## [1] 0.001Lower 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.
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.09849074902
##
## $f1
## [,1]
## [1,] 80.4430522
##
## $Q
## [,1]
## [1,] 4.724412174e-06
##
## $f
## [,1]
## [1,] 117.9467266
##
## $b
## [,1]
## [1,] 0.9217553011
##
## $mu0
## [1] 3.87576559e-05
##
## $theta
## [1] 0.1467324159
##
## $status
## [1] 5
##
## $LogLik
## [1] -623.8847506
##
## $objective
## [1] 623.8800548
##
## $message
## [1] "NLOPT_MAXEVAL_REACHED: Optimization stopped because maxeval (above) was reached."
##
## $limit
## [1] FALSE
##
## attr(,"class")
## [1] "spm.continuous"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.150267163523 0.001692834917
## [2,] 0.001979205983 -0.148012155129
##
## $f1
## [,1]
## [1,] 105.2777868
## [2,] 195.2059583
##
## $Q
## [,1] [,2]
## [1,] 1.302225923e-06 1.305236854e-07
## [2,] 1.299338548e-07 1.285194449e-06
##
## $f
## [,1]
## [1,] 107.4474372
## [2,] 210.4295117
##
## $b
## [,1]
## [1,] 1.119337586
## [2,] 1.958657445
##
## $mu0
## [1] 0.0001124533236
##
## $theta
## [1] 0.073700401
##
## $status
## [1] 5
##
## $LogLik
## [1] 1543.7251
##
## $objective
## [1] -1978.053682
##
## $message
## [1] "NLOPT_MAXEVAL_REACHED: Optimization stopped because maxeval (above) was reached."
##
## $limit
## [1] FALSE
##
## attr(,"class")
## [1] "spm.continuous"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## $a
## [1] -0.05
##
## $f1
## [1] 80
##
## $Q
## [1] 2e-08
##
## $f
## [1] 80
##
## $b
## [1] 5
##
## $mu0
## [1] 0.001Imagine 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## $a
## [1] -0.05
##
## $f1
## [1] 80
##
## $Q
## [1] 2e-08
##
## $b
## [1] 80
##
## $mu0
## [1] 5
##
## $<NA>
## <NA>
## NAAs you can see, there is no parameter \(f\) in \(opt.par\). This is 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)
# Temporarily commented below
# 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## $a
## [1] -0.05
##
## $f1
## [1] 80
##
## $Q
## [1] 2e-08
##
## $b
## [1] 5
##
## $mu0
## [1] 0.001You can do the same manner if you want two or more parameters to be equal to zero.
Function spm_con_1d(...) allows for very fast parameter estimating for one-dimensional model. This function implements a analytical solution to estimate the parameters in the continuous SPM model by assuming all the parameters are constants. Below there is an example.
library(stpm)
dat <- simdata_cont(N=500)
colnames(dat) <- c("id", "xi", "t1", "t2", "y", "y.next")
res <- spm_con_1d(as.data.frame(dat), a=-0.05, b=2, q=1e-8, f=80, f1=90, mu0=1e-3, theta=0.08)## [1] "Initial values:"
## [1] -5e-02 2e+00 1e-08 8e+01 9e+01 1e-03 8e-02
## [1] "Lower bounds:"
## [1] -2.5e-01 2.0e-01 1.0e-10 4.0e+01 4.5e+01 1.0e-04 8.0e-02
## [1] "Upper bounds:"
## [1] -2.5e-03 1.0e+01 1.0e-07 1.6e+02 1.8e+02 1.0e-02 8.0e-02res## $est
## Coeff. Std. Err. z p value
## a -0.05039551161 NaN NaN NaN
## b 4.98298558756 NaN NaN NaN
## q 0.00000000010 NaN NaN NaN
## f 79.99999939383 NaN NaN NaN
## f1 80.63419811679 NaN NaN NaN
## mu0 0.00010000000 NaN NaN NaN
## theta 0.08000000000 NaN NaN NaN
##
## $hessian
## a b q f
## a 2.683389364e+05 6.248489666e+03 NaN 6.371292106e-07
## b 6.248489666e+03 1.719845228e+03 NaN -1.179747026e-12
## q NaN NaN NaN NaN
## f 6.371292106e-07 -1.179747026e-12 NaN 5.805262619e-06
## f1 -2.377069205e+01 -2.230640043e-08 NaN 4.179718091e-08
## mu0 7.519472680e-01 2.124368232e-02 NaN -3.354831501e-03
## theta 1.124218513e-02 3.021276719e-04 NaN -5.674401848e-05
## f1 mu0 theta
## a -2.377069205e+01 7.519472680e-01 1.124218513e-02
## b -2.230640043e-08 2.124368232e-02 3.021276719e-04
## q NaN NaN NaN
## f 4.179718091e-08 -3.354831501e-03 -5.674401848e-05
## f1 3.268203320e+00 2.867614181e-04 4.880043763e-06
## mu0 2.867614181e-04 -2.090909102e+08 1.864276359e+09
## theta 4.880043763e-06 1.864276359e+09 1.720542968e+07
##
## $lik
## [1] 50734.94692
##
## $con
## [1] 4
##
## $message
## [1] "NLOPT_XTOL_REACHED: Optimization stopped because xtol_rel or xtol_abs (above) was reached."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).
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 86.15622307 80.76966086
## 2 1 0 31 32 80.76966086 81.78895997
## 3 1 0 32 33 81.78895997 75.01422405
## 4 1 0 33 34 75.01422405 67.21544319
## 5 1 0 34 35 67.21544319 69.52583209
## 6 1 0 35 36 69.52583209 67.69272745The 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 30.83416090 32.06491080 81.62963334 72.00396582
## 2 0 0 32.06491080 33.93339324 72.00396582 75.24141447
## 3 0 0 33.93339324 35.87695218 75.24141447 68.88479271
## 4 0 0 35.87695218 37.05707884 68.88479271 68.03964045
## 5 0 0 37.05707884 38.71625430 68.03964045 71.15197180
## 6 0 0 38.71625430 40.65291268 71.15197180 65.84308133Stochastic 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.
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 49.09835424 50.10808925 0 80.11502416 77.91797206
## 2 0 0 50.10808925 51.17717465 0 77.91797206 83.57870782
## 3 0 0 51.17717465 52.24380817 0 83.57870782 74.96832058
## 4 0 0 52.24380817 53.22381880 0 74.96832058 67.48488598
## 5 0 0 53.22381880 54.23775739 0 67.48488598 57.70916873
## 6 0 0 54.23775739 55.26161700 0 57.70916873 60.57452685In 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 80.69013266 81.59680529 0 79.51788742 80.94716668
## 2 0 0 81.59680529 82.49777652 0 80.94716668 87.06222956
## 3 0 0 82.49777652 83.44990544 0 87.06222956 80.43281572
## 4 0 0 83.44990544 84.47192504 0 80.43281572 83.88261902
## 5 0 0 84.47192504 85.50392549 0 83.88261902 84.18701602
## 6 0 0 85.50392549 86.55172258 0 84.18701602 82.52480198#Parameters estimation:
pars <- spm_pobs(x=data)## Parameter QL achieved lower/upper bound.
## 2.75e-08pars## $aH
## [,1]
## [1,] -0.05474260033
##
## $aL
## [,1]
## [1,] -0.01089591815
##
## $f1H
## [,1]
## [1,] 65.79736719
##
## $f1L
## [,1]
## [1,] 87.92559711
##
## $QH
## [,1]
## [1,] 2.149314265e-08
##
## $QL
## [,1]
## [1,] 2.75e-08
##
## $fH
## [,1]
## [1,] 56.06473469
##
## $fL
## [,1]
## [1,] 72.08992069
##
## $bH
## [,1]
## [1,] 3.9846025
##
## $bL
## [,1]
## [1,] 4.96746206
##
## $mu0H
## [1] 8.79222864e-06
##
## $mu0L
## [1] 9.003086071e-06
##
## $thetaH
## [1] 0.07203219785
##
## $thetaL
## [1] 0.09000455722
##
## $p
## [1] 0.2685599973
##
## $limit
## [1] TRUE
##
## attr(,"class")
## [1] "pobs.spm"Here and represents parameters when \(Z\) = 1 (H) and 0 (L).
library(stpm)
data.genetic <- sim_pobs(N=5, mode='observed')
head(data.genetic)## id xi t1 t2 Z y1 y1.next
## 1 0 0 95.81419462 96.88293917 0 80.47733468 78.56694203
## 2 0 0 96.88293917 97.90719645 0 78.56694203 73.63145533
## 3 0 0 97.90719645 98.90721046 0 73.63145533 75.20806657
## 4 0 0 98.90721046 99.82324608 0 75.20806657 81.01598947
## 5 0 0 99.82324608 100.85129431 0 81.01598947 80.90168024
## 6 0 0 100.85129431 101.82788991 0 80.90168024 80.54295148data.nongenetic <- sim_pobs(N=10, mode='unobserved')
head(data.nongenetic)## id xi t1 t2 y1 y1.next
## 1 0 0 56.46516128 57.41999032 79.16666008 68.11428709
## 2 0 0 57.41999032 58.41460089 68.11428709 63.19305721
## 3 0 0 58.41460089 59.37608884 63.19305721 58.18889797
## 4 0 0 59.37608884 60.28358280 58.18889797 57.20918105
## 5 0 0 60.28358280 61.20323114 57.20918105 53.87106167
## 6 0 0 61.20323114 62.13048508 53.87106167 53.58402084#Parameters estimation:
pars <- spm_pobs(x=data.genetic, y = data.nongenetic, mode='combined')
pars## $aH
## [,1]
## [1,] -0.05423840127
##
## $aL
## [,1]
## [1,] -0.00890013468
##
## $f1H
## [,1]
## [1,] 63.44067636
##
## $f1L
## [,1]
## [1,] 73.24680668
##
## $QH
## [,1]
## [1,] 1.911854551e-08
##
## $QL
## [,1]
## [1,] 2.466991874e-08
##
## $fH
## [,1]
## [1,] 64.03593701
##
## $fL
## [,1]
## [1,] 82.9550889
##
## $bH
## [,1]
## [1,] 4.398289506
##
## $bL
## [,1]
## [1,] 5.303978886
##
## $mu0H
## [1] 7.675561423e-06
##
## $mu0L
## [1] 9.03975955e-06
##
## $thetaH
## [1] 0.07242570039
##
## $thetaL
## [1] 0.09002392201
##
## $p
## [1] 0.265954788
##
## $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.
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.114382051head(ex_data$gene_data)## id geno
## 1 1 1
## 2 2 1
## 3 3 0
## 4 4 0
## 5 5 1
## 6 6 0res <- 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-07),
effect=c('q'), method = "tnewton")## [1] "Initial values:"
## [1] -2e-02 -2e-02 2e-01 2e-01 1e-02 1e-02 3e+00 3e+00 3e+00 3e+00
## [11] 1e-02 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+00
## [11] 1e-03 1e-03 1e-07
## [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-01 1e-05res## $est
## Coeff. Std. Err. z p value
## a -0.031238526396 9.063948538e-04 -3.446458932e+01 0.0000000000
## b 0.101329776535 2.785757812e-04 3.637422323e+02 0.0000000000
## q_0 0.000010000000 3.124724202e-03 3.200282442e-03 0.9974465484
## q_2 0.003608789876 5.278940497e-03 6.836201087e-01 0.4942150834
## f 2.999731527110 3.052250225e-02 9.827934494e+01 0.0000000000
## f1 3.000056947322 1.716252597e-02 1.748027622e+02 0.0000000000
## mu0 0.001000000000 1.552500478e+04 6.441221851e-08 0.9999999486
## theta 0.000000100000 8.987112435e+01 1.112704450e-09 0.9999999991
##
## $lik
## [1] -121546.402
##
## $con
## [1] -1
##
## $message
## [1] "NLOPT_FAILURE: Generic failure code."
##
## $hessian
## a b q_0 q_2
## a 2.354652039e+06 6.786921795e+05 158.058062204 82.3853495314
## b 6.786921795e+05 1.325446527e+07 1517.123161715 775.1749357615
## q_0 1.580580622e+02 1.517123162e+03 -39.193641214 -13.8346645758
## q_2 8.238534953e+01 7.751749358e+02 -13.834664576 -13.8346698668
## f 3.489487410e-01 -4.196232148e-03 36528.531980456 19000.8175924782
## f1 -8.446522178e+04 5.065699701e+01 -2.840122726 -1.4581382745
## mu0 2.228043741e-03 8.307251846e-04 1.053761684 0.6444883984
## theta 4.230830818e-01 1.744572073e-01 194.433867000 117.7827687934
## f f1 mu0 theta
## a 3.489487410e-01 -8.446522178e+04 2.228043741e-03 0.423083081841
## b -4.196232148e-03 5.065699701e+01 8.307251846e-04 0.174457207322
## q_0 3.652853198e+04 -2.840122726e+00 1.053761684e+00 194.433867000043
## q_2 1.900081759e+04 -1.458138275e+00 6.444883984e-01 117.782768793404
## f 1.670661873e+02 -1.202720644e-02 3.587745596e-03 0.465787947178
## f1 -1.202720644e-02 6.470562259e+03 -5.237961886e-05 -0.006925314665
## mu0 3.587745596e-03 -5.237961886e-05 0.000000000e+00 0.000000000000
## theta 4.657879472e-01 -6.925314665e-03 0.000000000e+00 0.000000000000
##
## $beta
## Coeff. Std. Err. Chi. Sq p value
## beta_a NA NA NA NA
## beta_b NA NA NA NA
## beta_q 0.001799394938 0.004201832138 0.0007705738917 0.9778541969
## beta_f NA NA NA NA
## beta_f1 NA NA NA NA
## beta_mu0 NA NA NA NAHere: 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.
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 NA
## 2 1 0 32 79.80944762
## 3 1 0 34 83.30071095
## 4 1 0 36 87.67442606
## 5 1 0 38 91.72665099
## 6 1 0 40 93.88187750##### Look at the imputed data #####
head(imp.data)## id xi t y1
## 1 1 0 30 91.88742030
## 2 1 0 32 79.80944762
## 3 1 0 34 83.30071095
## 4 1 0 36 87.67442606
## 5 1 0 38 91.72665099
## 6 1 0 40 93.88187750#########################################################
################ 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 NA NA
## 2 1 0 32 NA 23.61127799
## 3 1 0 34 79.12870268 26.60658498
## 4 1 0 36 NA 24.24929082
## 5 1 0 38 84.87254156 NA
## 6 1 0 40 89.53286827 32.41289417###### Look at the imputed data #####
head(imp.data)## id xi t y1 y2
## 1 1 0 30 95.46202972 25.66067721
## 2 1 0 32 95.23194281 23.61127799
## 3 1 0 34 79.12870268 26.60658498
## 4 1 0 36 79.79518739 24.24929082
## 5 1 0 38 84.87254156 23.87629038
## 6 1 0 40 89.53286827 32.41289417We provide a simple function to predict the next value of . Refer to the example below:
#library(stpm)
#data <- simdata_discr(N=100, format="long")
#res <- spm_discrete(data)
#splitted <- split(data, data$id)
#df <- data.frame()
#lapply(1:100, function(i) {df<<-rbind(df,splitted[[i]][dim(splitted[[i]])[1],c("id", "xi", "t1", "y1")])})
#names(df) <- c("id", "xi", "t", "y")
#predicted <- predict(object=res, data=df, dt=3)
#head(predicted)The package offers following five hypotheses to test for function (Arbeev et al. 2016):
H01: \(Q(t)=0\) (i.e., \(a_Q = 0\) and \(b_Q = 0\),so that there is no quadratic term in the hazard rate and mortality is described by the baseline Gompertz rate \(μ_0(t)\)).
H02: \(Q(t) = a_Q\) (i.e., \(b_Q = 0\)).
H03: \(f_1(t) = 0\) (i.e., \(a_{f1} = 0\) and \(b_{f1} = 0\)).
H04: \(f_1(t) = a_{f_1}\) (i.e., \(b_{f_1} = 0\)).
H05: \(a(t) = a_Y\) (i.e., \(b_Y = 0\)).
To perform hypothesis testing you should put the variable lrtest to TRUE (this is "H01" by default) or to any of the following: "H01", "H02", "H03", "H04", "H05".
library(stpm)
n <- 1000
# Data simulation:
data <- simdata_time_dep(N=n, format="long")
head(data)
# Hypotheses testing
## H01
res <- spm_time_dep(data, verbose=F,
frm = list(at="a", f1t="f1", Qt="Q", ft="f", bt="b", mu0t="mu0"),
start=list(a=-0.05, f1=80, Q=1e-8, f=90, b=5, mu0=0.001),
lb=c(a=-1, f1=30, Q=1e-9, f=10, b=1, mu0=1e-6),
ub=c(a=0, f1=120, Q=1e-7, f=150, b=10, mu0=1e-2),
opts = list(algorithm = "NLOPT_LN_NELDERMEAD",
maxeval = 200, ftol_rel = 1e-12), lrtest="H01")
res$alternative$lr.test.pval
## H02
res <- spm_time_dep(data, verbose=F,
frm = list(at="a", f1t="f1", Qt="1e-6", ft="f", bt="b", mu0t="mu0"),
start=list(a=-0.05, f1=80, Q=1e-8, f=90, b=5, mu0=0.001),
lb=c(a=-1, f1=30, Q=1e-9, f=10, b=1, mu0=1e-6),
ub=c(a=0, f1=120, Q=1e-7, f=150, b=10, mu0=1e-2),
opts = list(algorithm = "NLOPT_LN_NELDERMEAD",
maxeval = 200, ftol_rel = 1e-12), lrtest="H02")
res$alternative$lr.test.pval
## H03
res <- spm_time_dep(data, verbose=F,
frm = list(at="a", f1t="f1", Qt="Q", ft="f", bt="b", mu0t="mu0"),
start=list(a=-0.05, f1=80, Q=1e-8, f=90, b=5, mu0=0.001),
ub=c(a=0, f1=120, Q=1e-7, f=150, b=10, mu0=1e-2),
opts = list(algorithm = "NLOPT_LN_NELDERMEAD",
maxeval = 200, ftol_rel = 1e-12), lrtest="H03")
res$alternative$lr.test.pval
## H04
res <- spm_time_dep(data, verbose=F,
frm = list(at="a", f1t="120", Qt="Q", ft="f", bt="b", mu0t="mu0"),
start=list(a=-0.05, f1=80, Q=1e-8, f=90, b=5, mu0=0.001),
lb=list(a=-1, f1=30, Q=1e-9, f=10, b=1, mu0=1e-6),
opts = list(algorithm = "NLOPT_LN_NELDERMEAD",
maxeval = 200, ftol_rel = 1e-12), lrtest="H04")
res$alternative$lr.test.pval
## H05
res <- spm_time_dep(data, verbose=F,
frm = list(at="-0.1", f1t="f1", Qt="Q", ft="f", bt="b", mu0t="mu0"),
start=list(a=-0.05, f1=80, Q=1e-8, f=90, b=5, mu0=0.001),
opts = list(algorithm = "NLOPT_LN_NELDERMEAD",
maxeval = 200, ftol_rel = 1e-12), lrtest="H05")
res$alternative$lr.test.pvalAkushevich, 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.
Arbeev, Konstantin G., Alan A. Cohen, Liubov S. Arbeeva, Emmanuel Milot, Eric Stallard, Alexander M. Kulminski, Igor Akushevich, Svetlana V. Ukraintseva, Kaare Christensen, and Anatoliy I. Yashin. 2016. “Optimal Versus Realized Trajectories of Physiological Dysregulation in Aging and Their Relation to Sex-Specific Mortality Risk.” Frontiers in Public Health 4 (January). Frontiers Media SA. doi:10.3389/fpubh.2016.00003.
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.