inst/doc/stpm-vignette.R

In stpm: Stochastic Process Model for Analysis of Longitudinal and Time-to-Event Outcomes

## ---- message=FALSE, echo=FALSE, eval=FALSE-----------------------------------
#  library(knitcitations)
#  cleanbib()
#  options("citation_format" = "pandoc")
#  r<-citep("10.1016/0040-5809(77)90005-3")
#  r<-citep("10.1016/j.mbs.2006.11.006")
#  r<-citep("10.1080/08898480590932296")
#  r<-citep("10.1007/s10522-006-9073-3")
#  r<-citep("10.1016/j.jtbi.2009.01.023")
#  r<-citep("10.3389/fpubh.2014.00228")
#  r<-citep("10.1002/gepi.22058")
#  r<-citep("10.3389/fpubh.2016.00003")
#  write.bibtex(file="references.bib")

## ----eval=FALSE---------------------------------------------------------------
#  install.packages("stpm")

## ----eval=FALSE---------------------------------------------------------------
#  require(devtools)
#  devtools::install_github("izhbannikov/stpm")

## ----results='hide',warning=FALSE,echo=FALSE,message=FALSE--------------------
library(stpm)
# Reading longitude data:
longdat <- read.csv(system.file("extdata","longdat.csv",package="stpm"))

## ----echo=FALSE---------------------------------------------------------------
head(longdat)

## ---- echo=FALSE, message=FALSE-----------------------------------------------
data <- simdata_discr(N=1000, format="short")

## ---- echo=FALSE--------------------------------------------------------------
head(data)

## ---- echo=FALSE, message=FALSE-----------------------------------------------
data <- simdata_discr(N=1000, format="long")

## ---- echo=FALSE--------------------------------------------------------------
head(data)

## -----------------------------------------------------------------------------
library(stpm)
data <- simdata_discr(N=10) # simulate data for 10 individuals, "long" format (default)
head(data)

## -----------------------------------------------------------------------------
library(stpm)
data <- simdata_cont(N=5, format="short") # simulate data for 5 individuals, "short" format
head(data)

## -----------------------------------------------------------------------------
library(stpm)
#Data simulation (200 individuals)
data <- simdata_discr(N=100)
#Estimation of parameters
pars <- spm_discrete(data)
pars

## -----------------------------------------------------------------------------
library(stpm)
#Simulate some data for 50 individuals
data <- simdata_cont(N=50)
head(data)
#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

## -----------------------------------------------------------------------------
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

## -----------------------------------------------------------------------------
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

## -----------------------------------------------------------------------------
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


## -----------------------------------------------------------------------------
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



## -----------------------------------------------------------------------------
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

## -----------------------------------------------------------------------------
library(stpm)
set.seed(12345678)
n <- 10
data <- simdata_time_dep(N=n, format = "long")

# 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),
                        opts=list(maxit=100), 
                        verbose=F)
opt.par


## -----------------------------------------------------------------------------
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)
res

## -----------------------------------------------------------------------------
library(stpm)
data <- simdata_discr(N=10)
head(data)

## ---- eval=T------------------------------------------------------------------
library(stpm)
data <- simdata_cont(N=10)
head(data)

## ---- eval=T------------------------------------------------------------------
library(stpm)
data <- sim_pobs(N=10)
head(data)

## ---- eval=T------------------------------------------------------------------
library(stpm)
#Generating data:
data <- sim_pobs(N=10)
head(data)
#Parameters estimation:
pars <- spm_pobs(x=data)
pars

## ---- eval=T------------------------------------------------------------------
library(stpm)
data.genetic <- sim_pobs(N=5, mode='observed')
head(data.genetic)
data.nongenetic <- sim_pobs(N=10, mode='unobserved')
head(data.nongenetic)
#Parameters estimation:
pars <- spm_pobs(x=data.genetic, y = data.nongenetic, mode='combined')
pars

## ---- eval=T------------------------------------------------------------------
library(stpm) 
data(ex_spmcon1dg)
head(ex_data$spm_data)
head(ex_data$gene_data)
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-07), 
                    effect=c('q'), method = "tnewton")
res

## ---- eval=TRUE, message=FALSE------------------------------------------------
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)

##### Look at the imputed data #####
head(imp.data)

#########################################################
################ 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)

###### Look at the imputed data #####
head(imp.data)

## -----------------------------------------------------------------------------
#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)

## ---- eval=FALSE--------------------------------------------------------------
#  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.pval