R/DLSSM.R
In jiakunj/DLSSM: Dynamic Logistic State Space Prediction Model
Defines functions
DLSSM.plot
DLSSM
DLSSM.valid
DLSSM.smooth
DLSSM.filter
DLSSM.predict
DLSSM.init
Batched
Documented in
DLSSM.filter
DLSSM.init
DLSSM.plot
DLSSM.predict
DLSSM.smooth
DLSSM.valid
#' Combine data into Batched data
#' @author Jiakun Jiang, Wei Yang, Wensheng Guo
#' @description The time domain of observation will first be standardized into [0,1]. Then [0,1] will be divided into S equally spaced intervals as described in Jiang et al.(2021, Biometrics).
#' Then those intervals slice the dataset to S batches of data.
#' @param formula An object of class "formula" (or one that can be coerced to that class): a symbolic description of response and covariates in the model.
#' @param data Dataset matrix containing the observations (one rows is a sample).
#' @param time The time variable in the dataset. The varying coefficient functions are assumed to be smooth functions of this variable.
#' @param S Number of batches
#' @export
#' @return
#' \tabular{ll}{
#' \code{batched} \tab List of batched data, the element of list is matrix with each row representing a sample \cr
#' \tab \cr
#' \code{gap.len} \tab interval length 1/S \cr
#' }
Batched<-function(formula,data,time,S){
if(any(is.na(data))==T) { stop("The dataset have missing data") }
var_names=all.vars(formula)
if(all(c(var_names,time)%in%colnames(data))==F) {stop("There exist undefined variables in formula")}
y=data[,var_names[1]]
x=data[,var_names[-1]]
t=data[,time]
n=dim(x)[1] #sample size
q=dim(x)[2] #dimension of covariates
if(n!=length(t)|n!=length(y)){
stop("The dimension of data does not match.")
}
if(as.integer(abs(S))!=S){
stop("S need to be a positive integer")
}
t=t/max(t)
oder=order(t)
t.order=t[oder]
x.order=x[oder,]
y.order=y[oder]
# create batch
t.batch=list() #Observational time-points
x.batch=list() #Covariates corresponding to varying coefficients
y.batch=list() #Binary outcome
data.all=list()
Data.check=rep(NA,S)
for(i in 1:S){
sel=t.order>(i-1)/S & t.order<=i/S
t.batch[[i]]=t.order[sel]
x.batch[[i]]=x.order[sel,]
y.batch[[i]]=as.numeric(y.order[sel])
Data.check[i]=sum(sel)
mmm=cbind(y.batch[[i]],x.batch[[i]],t.batch[[i]])
colnames(mmm)=c(var_names,time)
data.all[[i]]=mmm
}
if(any(Data.check==0)){
stop("There exist at least one interval with no data")
}
invisible(list(batched=data.all,x.batch=x.batch,y.batch=y.batch,t.batch=t.batch,gap.len=1/S,S=S,time=time,formula=formula))
#invisible(list(x.batch=x.batch,y.batch=y.batch,t.batch=t.batch,gap.len=1/S))
}
#' Initial model fitting
#' @author Jiakun Jiang, Wei Yang and Wensheng Guo
#' @description This function is for tuning smoothing parameters using training data. The likelihood was calculated by Kalman Filter and maximized to estimate the smoothing parameters.
#' For the given smoothing parameters, the model coefficients can
#' be efficiently estimated using a Kalman filtering algorithm.
#' @param data.batched A object generated by function Data.batched()
#' @param S0 Number of batches of data to be used as training dataset
#' @param vary.effects The names of variables in the dataset assumed to have a time-varying regression effect on the outcome.
#' @param autotune T/F indicates whether or not the automatic tuning procedure described in Jiang et al. (2021) should be applied. Default is true.
#' @param Lambda Specify smoothing parameters if autotune=F
#' @import Matrix
#' @export
#' @importFrom "graphics" "abline" "lines" "par"
#' @importFrom "stats" "na.omit" "optim" "optimize"
#' @return
#' \tabular{ll}{
#' \code{Lambda:} \tab smoothing parameters \cr
#' \tab \cr
#' \code{Smooth:} \tab smoothed state vector \cr
#' \tab \cr
#' \code{Smooth.var:} \tab covariance of smoothed state vector in Smooth. \cr
#' }
#' @examples
#' \donttest{
#' set.seed(321)
#' n=8000
#' beta0=function(t) 0.1*t-1
#' beta1=function(t) cos(2*pi*t)
#' beta2=function(t) sin(2*pi*t)
#' alph1=alph2=1
#' x=matrix(runif(n*4,min=-4,max=4),nrow=n,ncol=4)
#' t=sort(runif(n))
#' coef=cbind(beta0(t),beta1(t),beta2(t),rep(alph1,n),rep(alph2,n))
#' covar=cbind(rep(1,n),x)
#' linear=apply(coef*covar,1,sum)
#' prob=exp(linear)/(1+exp(linear))
#' y=as.numeric(runif(n)<prob)
#' sim.data=cbind(y,x,t)
#' colnames(sim.data)=c("y","x1","x2","x3","x4","t")
#' formula = y~x1+x2+x3+x4
#' # Divide the time domain [0,1] into S=100 equally spaced intervals
#' S=100
#' S0=75
#' data.batched=Batched(formula, data=sim.data, time="t", S)
#'
#' # using first 75 batches as training dataset to tune smoothing parameters
#' fit0=DLSSM.init(data.batched, S0, vary.effects=c("x1","x2"))
#' fit0$Lambda
#' DLSSM.plot(fit0)
#' }
DLSSM.init<-function(data.batched,S0,vary.effects,autotune=TRUE,Lambda=NULL){
formula=data.batched$formula
S=data.batched$S
if(S<S0){stop("Number of batches of training dataset exceed S")}
x.names=all.vars(formula)[-1]
vary=match(vary.effects,x.names)
time=data.batched$time
x.b=data.batched$x.batch
y.b=data.batched$y.batch
t.b=data.batched$t.batch
q=dim(x.b[[1]])[2] #dimension of covariates filter
q1=num.vary=length(vary.effects) #number of covariates with varying coefficient
q2=q-q1 #number of covariates with time-invariant coefficient
if(num.vary>0&num.vary<length(x.names)){
t=list() # observational time-points
X=list() # Covariates corresponding to varying coefficients
XX=list() # Covariates corresponding to constant coefficients
y=list() # binary outcome
for(i in 1:S0){
t[[i]]=t.b[[i]]
X[[i]]=x.b[[i]][,vary]
XX[[i]]=x.b[[i]][,setdiff(1:q,vary)]
y[[i]]=as.numeric(y.b[[i]])
}
}
if(num.vary==0){ # no varying-coefficient
t=list()
XX=list()
y=list()
for(i in 1:S0){
t[[i]]=t.b[[i]]
XX[[i]]=x.b[[i]]
y[[i]]=as.numeric(y.b[[i]])
}
}
if(num.vary==length(x.names)){
t=list() # observational time-points
X=list() # Covariates corresponding to varying coefficients
y=list() # binary outcome
for(i in 1:S0){
t[[i]]=t.b[[i]]
X[[i]]=x.b[[i]][,vary]
y[[i]]=as.numeric(y.b[[i]])
}
}
dim=q1+1 # Number of varying coefficients including intercept
dim.con=q2 # Number of constant coefficients
TT=matrix(0,2*dim+dim.con,2*dim+dim.con) # Transformation matrix in state-space model
QQ=matrix(NA,2,2) # Covariance matrix in state-space model
delta.t=data.batched$gap.len # Batch data be generated equally spaced
TTq1=matrix(0,2*dim,2*dim) # Transformation matrix for varying coefficient
TTq2=diag(rep(1,dim.con)) # Transformation matrix for constant coefficient
for(ss in 1:dim){
TTq1[(2*ss-1):(2*ss),(2*ss-1):(2*ss)]=matrix(c(1, 0, delta.t, 1),2, 2)
}
TT=as.matrix(bdiag(TTq1,TTq2))
QQ=matrix(c((delta.t)^3/3,(delta.t)^2/2,(delta.t)^2/2,delta.t),2,2)
a1=rep(0,2*dim+dim.con) # Initial value
diff=log(100) # Diffuse variance of initial prior distribution
n.train=S0 # number of batches S
#########################################
#Likelihood function to tune smoothing parameters Lambda based on training data
#########################################
if(autotune==TRUE){
likehood.predictive=function(xx) {
Lambda=xx
Q=matrix(0,2*dim+dim.con,2*dim+dim.con)
Q1=matrix(0,2*dim,2*dim)
for(ss in 1:dim){
Q1[(2*ss-1):(2*ss),(2*ss-1):(2*ss)]=exp(Lambda[ss])*QQ
}
Q[1:(2*dim),1:(2*dim)]=Q1
P1 <- exp(diff)*diag(2*dim+dim.con)
Sample.likehood=rep(NA,n.train)
prediction=matrix(NA,n.train,dim*2+dim.con)
prediction.var=array(NA,dim=c(n.train,dim*2+dim.con,dim*2+dim.con))
prediction[1,]=a1
prediction.var[1,,]=P1
filter=matrix(NA,n.train,dim*2+dim.con)
filter.var=array(NA,dim=c(n.train,dim*2+dim.con,dim*2+dim.con))
intial=prediction[1,]
ZZ=matrix(0,2*dim+dim.con,length(t[[1]]))
if(dim>1&dim<(length(x.names)+1)){
ZZ[1,]=1
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[1]])[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[1]])
}
if(dim==1){
ZZ[1,]=1
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[1]])
}
if(dim==(length(x.names)+1)){
ZZ[1,]=1
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[1]])[sss,]
}
#ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[1]])
}
for(inter in 1:1000){
summa=rep(0,2*dim+dim.con)
summa.cov=matrix(0,2*dim+dim.con,2*dim+dim.con)
for(j in 1:length(t[[1]])){
y.hat=exp(ZZ[,j]%*%intial)/(1+exp(ZZ[,j]%*%intial))
summa=summa+as.numeric(y[[1]][j]-y.hat)*ZZ[,j]
summa.cov=summa.cov+as.numeric(y.hat*(1-y.hat))*ZZ[,j]%o%ZZ[,j]
}
SCORE=summa-t(solve(prediction.var[1,,])%*%(intial-matrix(prediction[1,],2*dim+dim.con,1)))
COVMATRIX=-solve(prediction.var[1,,])-summa.cov
recur=intial-(0.2*solve(COVMATRIX)%*%t(SCORE))
befor=intial
intial=recur
if(mean(abs(befor-recur))<0.00001){break}
}
filter[1,]=recur
filter.var[1,,]=solve(-COVMATRIX)
lihood=0
for(j in 1:length(t[[1]])){
y.hatt=exp(ZZ[,j]%*%filter[1,])/(1+exp(ZZ[,j]%*%filter[1,]))
Covmatrix=-solve(prediction.var[1,,])-as.numeric(y.hatt*(1-y.hatt))*ZZ[,j]%o%ZZ[,j]
densi=1/sqrt(det(prediction.var[1,,]))*
exp(-0.5*(filter[1,]-prediction[1,])%*%solve(prediction.var[1,,])%*%(filter[1,]-prediction[1,]))
lihood=lihood+log(2*pi*sqrt(det(solve(-Covmatrix)))*(y.hatt^(y[[1]][j]))*((1-y.hatt)^(1-y[[1]][j]))*densi)
}
Sample.likehood[1]=lihood
for(i in 2:n.train) {
prediction[i,]=TT%*%filter[i-1,]
prediction.var[i,,]=TT%*%filter.var[i-1,,]%*%t(TT)+Q
intial=prediction[i,]
ZZ=matrix(0,dim*2+dim.con,length(t[[i]]))
ZZ[1,]=1 #num.vary>0&num.vary<length(x.names)
if(num.vary>0&num.vary<length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
if(num.vary==0){
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
if(num.vary==length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
}
for(inter in 1:1000){
summa=rep(0,dim*2+dim.con)
summa.cov=matrix(0,dim*2+dim.con,dim*2+dim.con)
for(j in 1:length(t[[i]])){
y.hat=exp(ZZ[,j]%*%intial)/(1+exp(ZZ[,j]%*%intial))
summa=summa+as.numeric(y[[i]][j]-y.hat)*ZZ[,j]
summa.cov=summa.cov+as.numeric(y.hat*(1-y.hat))*ZZ[,j]%o%ZZ[,j]
}
SCORE=summa-t(solve(prediction.var[i,,])%*%(intial-matrix(prediction[i,],2*dim+dim.con,1)))
COVMATRIX=-solve(prediction.var[i,,])-summa.cov
recur=intial-(0.2*solve(COVMATRIX)%*%t(SCORE))
befor=intial
intial=recur
if(max(abs(befor-recur))<0.00001){break}
}
filter[i,]=recur
filter.var[i,,]=solve(-COVMATRIX)
lihood=0
for(j in 1:length(t[[i]])){
y.hatt=exp(ZZ[,j]%*%filter[i,])/(1+exp(ZZ[,j]%*%filter[i,]))
Covmatrix=-solve(prediction.var[i,,])-as.numeric(y.hatt*(1-y.hatt))*ZZ[,j]%o%ZZ[,j]
densi=1/sqrt(det(prediction.var[i,,]))*
exp(-0.5*(filter[i,]-prediction[i,])%*%solve(prediction.var[i,,])%*%(filter[i,]-prediction[i,]))
lihood=lihood+log(2*pi*sqrt(det(solve(-Covmatrix)))*(y.hatt^(y[[i]][j]))*((1-y.hatt)^(1-y[[i]][j]))*densi)
}
Sample.likehood[i]=lihood
}
likelihood=-sum((Sample.likehood))
return(likelihood)
}
if(dim>1){
search.opt=optim(rep(2,dim),likehood.predictive) #Optimize likelihood function with initial value (0,0)
Lambda=exp(search.opt$par)
}
if(dim==1){
search.opt=optimize(likehood.predictive,c(-20,20)) #Optimize likelihood function with initial value (0,0)
Lambda=exp(search.opt$minimum)
}
#Smoothing parameter for intercept
}else{
if(length(Lambda)!=num.vary+1){
stop("The number of smoothing parameters should equal to number of varying coefficients, including varying intercept.")
}
Lambda=Lambda
}
P1<-exp(diff)*diag(dim*2+dim.con) #Diffuse prior
Q=matrix(0,dim*2+dim.con,dim*2+dim.con)
Q1=matrix(0,2*dim,2*dim)
for(ss in 1:dim){
Q1[(2*ss-1):(2*ss),(2*ss-1):(2*ss)]=Lambda[ss]*QQ
}
Q[1:(2*dim),1:(2*dim)]=Q1
#########################################
#Model estimation using all data with Kalman filter
#########################################
prediction=matrix(NA,S0,2*dim+dim.con) # Mean prediction
prediction.var=array(NA,dim=c(S0,2*dim+dim.con,2*dim+dim.con)) # Covariance prediction
prediction[1,]=a1 # Initial mean
prediction.var[1,,]=P1 # Initial covariance
filter=matrix(NA,S0,2*dim+dim.con) # Filtering mean
filter.var=array(NA,dim=c(S0,2*dim+dim.con,2*dim+dim.con)) # Filtering covariance
#Kalman filter on first timepoint
intial=prediction[1,]
ZZ=matrix(0,2*dim+dim.con,length(t[[1]])) #define covariate z in our state space model (5)
ZZ[1,]=1
if(num.vary>0&num.vary<length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[1]])[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[1]])
}
if(num.vary==0){
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[1]])
}
if(num.vary==length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[1]])[sss,]
}
}
for(inter in 1:1000){ # Newton-Raphson iterative algorithm to get filtering estimator
summa=rep(0,2*dim+dim.con)
summa.cov=matrix(0,2*dim+dim.con,2*dim+dim.con)
for(j in 1:length(t[[1]])){
y.hat=exp(ZZ[,j]%*%intial)/(1+exp(ZZ[,j]%*%intial))
summa=summa+as.numeric(y[[1]][j]-y.hat)*ZZ[,j]
summa.cov=summa.cov+as.numeric(y.hat*(1-y.hat))*ZZ[,j]%o%ZZ[,j]
}
SCORE=summa-t(solve(prediction.var[1,,])%*%(intial-matrix(prediction[1,],2*dim+dim.con,1)))
COVMATRIX=-solve(prediction.var[1,,])-summa.cov
recur=intial-(0.2*solve(COVMATRIX)%*%t(SCORE))
befor=intial
intial=recur
if(mean(abs(befor-recur))<0.00001){break}
if(inter==1000){stop("failed converge","\n")}
}
filter[1,]=recur # Filtering mean
filter.var[1,,]=solve(-COVMATRIX) # Filtering covariance
#Kalman filter on timepoints from 2 to S
Prob.est=list()
for(i in 2:S0) {
# Prediction step
prediction[i,]=TT%*%filter[i-1,]
prediction.var[i,,]=TT%*%filter.var[i-1,,]%*%t(TT)+Q
# Filtering step
intial=prediction[i,]
ZZ=matrix(0,dim*2+dim.con,length(t[[i]]))
ZZ[1,]=1
if(num.vary>0&num.vary<length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
if(num.vary==0){
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
if(num.vary==length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
}
if(i>n.train){
y.hat=exp(t(ZZ)%*%intial)/(1+exp(t(ZZ)%*%intial))
Prob.est[[i-n.train]]=y.hat
}
for(inter in 1:1000){ # Newton-Raphson
summa=rep(0,dim*2+dim.con)
summa.cov=matrix(0,dim*2+dim.con,dim*2+dim.con)
for(j in 1:length(t[[i]])){
y.hat=exp(ZZ[,j]%*%intial)/(1+exp(ZZ[,j]%*%intial))
summa=summa+as.numeric(y[[i]][j]-y.hat)*ZZ[,j]
summa.cov=summa.cov+as.numeric(y.hat*(1-y.hat))*ZZ[,j]%o%ZZ[,j]
}
SCORE=summa-t(solve(prediction.var[i,,])%*%(intial-matrix(prediction[i,],2*dim+dim.con,1)))
COVMATRIX=-solve(prediction.var[i,,])-summa.cov
recur=intial-(0.2*solve(COVMATRIX)%*%t(SCORE))
befor=intial
intial=recur
if(mean(abs(befor-recur))<0.00001){break}
if(inter==1000){stop("Newton-Raphson failed to converge","\n")}
#if(inter==1000){break}
}
filter[i,]=recur
filter.var[i,,]=solve(-COVMATRIX)
}
#Smoothed coefficients on training data
state.smooth=matrix(NA,S0,2*dim+dim.con)
F.t.smooth=array(NA,dim=c(S0,2*dim+dim.con,2*dim+dim.con))
state.smooth[S0,]=filter[S0,]
F.t.smooth[S0,,]=filter.var[S0,,]
for(i in (S0-1):1){
A.i=filter.var[i,,]%*%t(TT)%*%solve(prediction.var[i+1,,])
state.smooth[i,]=filter[i,]+A.i%*%(state.smooth[i+1,]-prediction[i+1,])
F.t.smooth[i,,]=filter.var[i,,]-A.i%*%(prediction.var[i+1,,]-F.t.smooth[i+1,,])%*%t(A.i)
}
Est=list(pred=prediction,pred.var=prediction.var,filter=filter
,filter.var=filter.var,smooth=state.smooth,smooth.var=F.t.smooth
,Lambda=Lambda,vary.effects=vary.effects,q=q,q1=q1,q2=q2
,dim=dim,dim.con=dim.con,num.vary=num.vary,TT=TT,Q=Q,time=time
,Prob.est=as.vector(Prob.est),S=S,S0=S0,gap.len=delta.t,formula=formula,initial.fit=TRUE)
return(Est)
}
DLSSM.predict<-function(fit,K,newx){
vary=fit$vary
num.vary=length(vary)
dim=fit$dim
dim.con=fit$dim.con
if(is.vector(fit$filter)==T){
fil.last=fit$filter
filt.var.last=fit$filter.var
}
else{
fil.last=fit$filter[dim(fit$filter)[1],]
filt.var.last=fit$filter.var[dim(fit$filter)[1],,]
}
coef.prediction=matrix(NA,K,2*dim+dim.con) # Mean prediction
coef.prediction.var=array(NA,dim=c(K,2*dim+dim.con,2*dim+dim.con)) # Covariance prediction
for(pred in 1:K){
TT=fit$TT
Q=fit$Q
coef.prediction[pred,]=TT%*%fil.last
coef.prediction.var[pred,,]=TT%*%filt.var.last%*%t(TT)+Q
}
coef.pred.onlyK=coef.prediction[K,] # only k-step
coef.pred.onlyK.var=coef.prediction.var[K,,]
if(is.null(newx)==TRUE){prob.pred=NULL}
if(is.null(newx)==FALSE){
if(is.vector(newx)==T) {newx=matrix(newx,nrow=1)}
prediction=coef.pred.onlyK
x.names=all.vars(fit$formula)[-1]
if(all(c(x.names)%in%colnames(newx))==F) {stop("There exist variable in model but not in data")}
newx=newx[,x.names]
vary.effects=fit$vary.effects
vary=match(vary.effects,x.names[-1])
q=dim(newx)[2]
newx=newx[,c(vary,setdiff(1:q,vary))]
ZZ=matrix(0,dim(newx)[1],dim*2+dim.con)
ZZ[,1]=1
for(sss in 1:(dim-1)){
ZZ[,2*sss+1]=newx[,sss]
}
ZZ[,(2*dim+1):(2*dim+dim.con)]=newx[,dim:fit$q]
prob.pred=exp(ZZ%*%prediction)/(1+exp(ZZ%*%(prediction)))
}
return(list(coef.pred=coef.pred.onlyK,coef.pred.var=coef.pred.onlyK.var,prob.pred=prob.pred))
}
DLSSM.filter<-function(fit,newdata){
formula=fit$formula
x.names=all.vars(formula)[-1]
time=fit$time
if(all(c(x.names,time)%in%colnames(newdata))==F) {stop("There exist variable in model but not in data")}
newx=newdata[,x.names]
newy=newdata[,all.vars(formula)[1]]
vary.effects=fit$vary.effects
vary=match(vary.effects,x.names)
q=dim(newx)[2]
newx=newx[,c(vary,setdiff(1:q,vary))]
num.vary=length(vary.effects)
dim=fit$dim
dim.con=fit$dim.con
TT=fit$TT
Q=fit$Q
if(is.vector(fit$filter)==T){
fil.last=fit$filter
filt.var.last=fit$filter.var
}
if(is.vector(fit$filter)==F){
fil.last=fit$filter[dim(fit$filter)[1],]
filt.var.last=fit$filter.var[dim(fit$filter)[1],,]
}
predict=TT%*%fil.last
predict.var=TT%*%filt.var.last%*%t(TT)+Q
# Filtering step
intial=predict
ZZ=matrix(0,2*dim+dim.con,length(newy))
ZZ[1,]=1#num.vary>0&num.vary<length(x.names)
if(num.vary>0&num.vary<length(x.names)){
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(newx)[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(newx)[dim:q,]
}
if(num.vary==0){
ZZ[(2*dim+1):(2*dim+dim.con),]=t(newx)
}
if(num.vary==length(x.names)){
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(newx)[sss,]
}
#ZZ[(2*dim+1):(2*dim+dim.con),]=t(newx)[dim:q,]
}
for(inter in 1:1000){ # Newton-Raphson
summa=rep(0,dim*2+dim.con)
summa.cov=matrix(0,dim*2+dim.con,dim*2+dim.con)
for(j in 1:length(newy)){
y.hat=exp(ZZ[,j]%*%intial)/(1+exp(ZZ[,j]%*%intial))
summa=summa+as.numeric(newy[j]-y.hat)*ZZ[,j]
summa.cov=summa.cov+as.numeric(y.hat*(1-y.hat))*ZZ[,j]%o%ZZ[,j]
}
SCORE=summa-t(solve(predict.var)%*%(intial-matrix(predict,2*dim+dim.con,1)))
COVMATRIX=-solve(predict.var)-summa.cov
recur=intial-(0.2*solve(COVMATRIX)%*%t(SCORE))
befor=intial
intial=recur
if(mean(abs(befor-recur))<0.00001){break}
if(inter==1000){stop("failed converge","\n")}
}
filter=as.vector(recur)
filter.var=solve(-COVMATRIX)
Lambda=fit$Lambda
q=fit$q
q1=fit$q1
Est=list(filter=filter,filter.var=filter.var,Lambda=Lambda,vary=vary,q=q,q1=q1,dim=dim,dim.con=dim.con,TT=TT,Q=Q,time=time,formula=formula,vary.effects=vary.effects)
return(Est)
}
DLSSM.smooth<-function(fit,prediction,prediction.var,filter,filter.var){
dim=fit$dim
dim.con=fit$dim.con
TT=fit$TT
n.train=dim(filter)[1]
state.smooth=matrix(NA,n.train,2*dim+dim.con)
F.t.smooth=array(NA,dim=c(n.train,2*dim+dim.con,2*dim+dim.con))
state.smooth[n.train,]=filter[n.train,]
F.t.smooth[n.train,,]=filter.var[n.train,,]
for(i in (n.train-1):1){
A.i=filter.var[i,,]%*%t(TT)%*%solve(prediction.var[i+1,,])
state.smooth[i,]=filter[i,]+A.i%*%(state.smooth[i+1,]-prediction[i+1,])
F.t.smooth[i,,]=filter.var[i,,]-A.i%*%(prediction.var[i+1,,]-F.t.smooth[i+1,,])%*%t(A.i)
}
return(list(smooth=state.smooth,smooth.var=F.t.smooth))
}
#' Dynamical prediction on validation dataset
#' @author Jiakun Jiang, Wei Yang and Wensheng Guo
#' @description After we have fitted initial model, we can do validation. It is iteratively doing K-steps ahead prediction and model updating (filtering)
#' when a new batch of data becomes available.
#' The validation include K-steps ahead prediction of state vector and probabilities on validation interval.
#' @param fit0 Initial fitted model
#' @param data.batched Batched dataset generated by function Batched()
#' @param K Number of steps for ahead prediction
#' @import Matrix
#' @export
#' @details The argument fit could be object of DLSSM or DLSSM.init.
#' @importFrom "graphics" "abline" "lines" "par"
#' @importFrom "stats" "na.omit" "optim" "optimize"
#' @return
#' \tabular{ll}{
#' \code{pred.K:} \tab K-steps ahead predicted coefficients \cr
#' \tab \cr
#' \code{pred.var.K:} \tab covariance of K-steps ahead predicted coefficients \cr
#' \tab \cr
#' \code{pred.prob.K:} \tab K-steps ahead predicted probabilities \cr
#' }
#' @examples
#'\donttest{
#' set.seed(321)
#' n=8000
#' beta0=function(t) 0.1*t-1
#' beta1=function(t) cos(2*pi*t)
#' beta2=function(t) sin(2*pi*t)
#' alph1=alph2=1
#' x=matrix(runif(n*4,min=-4,max=4),nrow=n,ncol=4)
#' t=sort(runif(n))
#' coef=cbind(beta0(t),beta1(t),beta2(t),rep(alph1,n),rep(alph2,n))
#' covar=cbind(rep(1,n),x)
#' linear=apply(coef*covar,1,sum)
#' prob=exp(linear)/(1+exp(linear))
#' y=as.numeric(runif(n)<prob)
#' sim.data=cbind(y,x,t)
#' colnames(sim.data)=c("y","x1","x2","x3","x4","t")
#' formula = y~x1+x2+x3+x4
#' # Divide the time domain [0,1] into S=100 equally spaced intervals
#' S=100
#' S0=75
#' data.batched=Batched(formula, data=sim.data, time="t", S)
#'
#' # using first 75 batches as training dataset to tune smoothing parameters
#' fit0=DLSSM.init(data.batched, S0, vary.effects=c("x1","x2"))
#' fit0$Lambda
#'
#' #After initial model fitting on training data, we move to dynamic prediction
#' fit=DLSSM.valid(fit0, data.batched, K=1)
#' DLSSM.plot(fit)
#' }
DLSSM.valid<-function(fit0,data.batched,K){
S=fit0$S
S0=fit0$S0
#if(S0+K>S){stop("S0+K>S happen, the predicted horizon exceed the whole dataset")}
formula=fit0$formula
vary.effects=fit0$vary.effects
TT=fit0$TT
Q=fit0$Q
x.names=all.vars(formula)[-1]
vary=match(vary.effects,x.names)
x.b=data.batched$x.batch
y.b=data.batched$y.batch
t.b=data.batched$t.batch
q=dim(x.b[[1]])[2] #dimension of covariates filter
q1=num.vary=length(vary.effects) #number of covariates with varying coefficient
q2=q-q1 #number of covariates with time-invariant coefficient
if(num.vary>0&num.vary<length(x.names)){
t=list() # observational time-points
X=list() # Covariates corresponding to varying coefficients
XX=list() # Covariates corresponding to constant coefficients
y=list() # binary outcome
for(i in 1:S){
t[[i]]=t.b[[i]]
X[[i]]=x.b[[i]][,vary]
XX[[i]]=x.b[[i]][,setdiff(1:q,vary)]
y[[i]]=as.numeric(y.b[[i]])
}
}
if(num.vary==0){ # no varying-coefficient covariates
t=list()
XX=list()
y=list()
for(i in 1:S){
t[[i]]=t.b[[i]]
XX[[i]]=x.b[[i]]
y[[i]]=as.numeric(y.b[[i]])
}
}
if(num.vary==length(x.names)){
t=list() # observational time-points
X=list() # Covariates corresponding to varying coefficients
#XX=list() # Covariates corresponding to constant coefficients
y=list() # binary outcome
for(i in 1:S){
t[[i]]=t.b[[i]]
X[[i]]=x.b[[i]][,vary]
#XX[[i]]=x.b[[i]][,setdiff(1:q,vary)]
y[[i]]=as.numeric(y.b[[i]])
}
}
dim=fit0$dim
dim.con=fit0$dim.con
dimen=2*dim+dim.con
Prediction=matrix(NA,S,dimen) # 3 varying-coefficients and 2 constants coefficients
Prediction.var=array(NA,dim=c(S,dimen,dimen))
Filtering=matrix(NA,S,dimen)
Filtering.var=array(NA,dim=c(S,dimen,dimen))
Prediction[1:S0,]=fit0$pred
Prediction.var[1:S0,,]=fit0$pred.var
Filtering[1:S0,]=fit0$filter
Filtering.var[1:S0,,]=fit0$filter.var
fit.last=fit0
for(i in (S0+1):S){
pred=DLSSM.predict(fit.last,newx=NULL,K=1)
Prediction[i,]=pred$coef.pred
Prediction.var[i,,]=pred$coef.pred.var
fit.last=DLSSM.filter(fit.last,data.batched$batched[[i]])
Filtering[i,]=fit.last$filter
Filtering.var[i,,]=fit.last$filter.var
}
prediction.K=matrix(NA,S,2*dim+dim.con) #for K steps ahead coef prediction
prediction.var.K=array(NA,dim=c(S,2*dim+dim.con,2*dim+dim.con))
Prob.est.K=list() # K steps ahead prediction
for(i in (S0+K):S){
K.S.P=matrix(NA,K,2*dim+dim.con)
K.S.P.var=array(NA,dim=c(K,2*dim+dim.con,2*dim+dim.con))
K.S.P[1,]=TT%*%Filtering[i-K,]
K.S.P.var[1,,]=TT%*%Filtering.var[i-K,,]%*%t(TT)+Q
if(K>1){
for(kp in 2:K){
K.S.P[kp,]=TT%*%K.S.P[kp-1,]
K.S.P.var[kp,,]=TT%*%K.S.P.var[kp-1,,]%*%t(TT)+Q
}
}
prediction.K[i,]=K.S.P[K,]
prediction.var.K[i,,]=K.S.P.var[K,,]
#K-steps ahead prediction of probability.
y.hat.k=list()
intial=prediction.K[i,]
ZZ=matrix(0,dim*2+dim.con,length(t[[i]]))
if(num.vary==q){
ZZ[1,]=1
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
}
if(num.vary>0&num.vary<q){
ZZ[1,]=1
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
if(num.vary==0){
ZZ[1,]=1
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
y.hat.k=exp(t(ZZ)%*%intial)/(1+exp(t(ZZ)%*%intial))
Prob.est.K[[i]]=y.hat.k
}
Smoothed=DLSSM.smooth(fit.last,Prediction,Prediction.var,Filtering,Filtering.var)
Est=list(pred=Prediction,pred.var=Prediction.var,filter=Filtering
,filter.var=Filtering.var,smooth=Smoothed$smooth,smooth.var=Smoothed$smooth.var
,pred.K=prediction.K[(S0+1):S,],pred.var.K=prediction.var.K[(S0+1):S,,],pred.prob.K=as.vector(Prob.est.K),Lambda=fit0$Lambda,vary.effects=vary.effects,q=q,q1=q1,q2=q2
,dim=dim,dim.con=dim.con,TT=TT,Q=Q
,S=S,S0=S0,gap.len=fit0$gap.len,K=K,formula=formula,initial.fit=FALSE,fit0=fit0)
return(Est)
}
#' Combine model training and validation in a integrated function
#' @author Jiakun Jiang, Wei Yang and Wensheng Guo
#' @description This combine model training and validation in a integrated automatic function DLSSM().
#' @param data.batched A object generated by function Data.batched()
#' @param S0 Number of batches of data to be used as training dataset
#' @param vary.effects The names of variables in the dataset assumed to have a time-varying regression effect on the outcome.
#' @param autotune T/F indicates whether or not the automatic tuning procedure desribed in Jiakun et al. (2021) should be applied. Default is true.
#' @param Lambda Specify smoothing parameters if autotune=F
#' @param K Number of steps for ahead prediction
#' @import Matrix
#' @importFrom "graphics" "abline" "lines" "par"
#' @importFrom "stats" "na.omit" "optim" "optimize"
#' @export
#' @return
#' \tabular{ll}{
#' \code{Lambda:} \tab smoothing parameters \cr
#' \tab \cr
#' \code{Smooth:} \tab smoothed state vector \cr
#' \tab \cr
#' \code{Smooth.var:} \tab covariance of smoothed state vector in Smooth. \cr
#' }
#' @examples
#' \donttest{
#' set.seed(321)
#' n=8000
#' beta0=function(t) 0.1*t-1
#' beta1=function(t) cos(2*pi*t)
#' beta2=function(t) sin(2*pi*t)
#' alph1=alph2=1
#' x=matrix(runif(n*4,min=-4,max=4),nrow=n,ncol=4)
#' t=sort(runif(n))
#' coef=cbind(beta0(t),beta1(t),beta2(t),rep(alph1,n),rep(alph2,n))
#' covar=cbind(rep(1,n),x)
#' linear=apply(coef*covar,1,sum)
#' prob=exp(linear)/(1+exp(linear))
#' y=as.numeric(runif(n)<prob)
#' sim.data=cbind(y,x,t)
#' colnames(sim.data)=c("y","x1","x2","x3","x4","t")
#' formula = y~x1+x2+x3+x4
#' # Divide the time domain [0,1] into S=100 equally spaced intervals
#' S=100
#' S0=75
#' data.batched=Batched(formula, data=sim.data, time="t", S)
#'
#'# Take first S0=75 batches as training data, remaining S-S0=25 batches of data as validation data.
#' fit1=DLSSM(data.batched, S0, vary.effects=c("x1","x2"), autotune=TRUE, Lambda=NULL, K=1)
#' DLSSM.plot(fit1)
#' fit2=DLSSM(data.batched, S0, vary.effects=c("x1","x2"), autotune=TRUE, Lambda=NULL, K=2)
#' DLSSM.plot(fit2)
#' }
DLSSM<-function(data.batched, S0, vary.effects, autotune=TRUE, Lambda=NULL, K){
fit0=DLSSM.init(data.batched,S0,vary.effects=vary.effects,autotune,Lambda)
fit=DLSSM.valid(fit0, data.batched, K)
Est=list(pred=fit0$pred,pred.var=fit0$pred.var,filter=fit0$filter
,filter.var=fit0$filter.var,smooth=fit0$smooth,smooth.var=fit0$smooth.var
,pred.K=fit$pred.K,pred.var.K=fit$pred.var.K,pred.prob.K=fit$pred.prob.K,Lambda=fit0$Lambda,vary.effects=vary.effects,q=fit$q,q1=fit$q1,q2=fit$q2
,dim=fit$dim,dim.con=fit$dim.con,TT=fit$TT,Q=fit$Q
,S=fit$S,S0=S0,gap.len=fit0$gap.len,K=K,formula=fit$formula,initial.fit=FALSE,fit0=fit0)
return(Est)
}
#' Plot coefficients
#' @author Jiakun Jiang, Wei Yang and Wensheng Guo
#' @description Plot smoothed coefficients in the training part and predicted coefficients in validation part, the two parts are divided by vertical dash line.
#' @param fit fitted object
#' @importFrom "graphics" "abline" "lines" "par"
#' @importFrom "stats" "na.omit" "optim" "optimize"
#' @export
#' @return Figures
#' @details If argument "fit" is an initial fitted model then only smoothed coefficients part are plotted.
DLSSM.plot<-function(fit){
x.names=all.vars(fit$formula)[-1]
vary.names=c(0,fit$vary.effects)
con.names=setdiff(x.names,vary.names)
S=fit$S
S0=fit$S0
q2=fit$q2
K=fit$K
initial.fit=fit$initial.fit
rows=length(fit$vary)+1
training_samp=fit$S0
cc=rows#+q2
c1=ceiling(sqrt(cc))
f1=floor(sqrt(cc))
if(cc==c1*f1){numb.plot=c(f1,c1)}
if(c1*f1>cc){numb.plot=c(f1,c1)}
if(c1*f1<cc){numb.plot=c(c1,c1)}
#if(rows>1&rows<length(x.names)+1){
plot.index=2*(1:rows)-1
#}
#if(rows==length(x.names)+1){
# plot.index=2*(1:rows)-1
#}
#if(rows==1){
# plot.index=2*(1:rows)-1
#}
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mfrow=numb.plot,mar=c(4, 4, 1, 1),oma=c(1,1,1,1))
scal=3
if(S0<S){
if(initial.fit==F){
est.p=fit$pred.K
est.var.p=fit$pred.var.K
fit0=fit$fit0
est.s=fit0$smooth
est.var.s=fit0$smooth.var
for(ss in plot.index){
p.up1=est.p[1:(S-S0-K+1),ss]+2*sqrt(est.var.p[1:(S-S0-K+1),ss,ss])
p.low1=est.p[1:(S-S0-K+1),ss]-2*sqrt(est.var.p[1:(S-S0-K+1),ss,ss])
s.up1=est.s[(1:S0),ss]+2*sqrt(est.var.s[(1:S0),ss,ss])
s.low1=est.s[(1:S0),ss]-2*sqrt(est.var.s[(1:S0),ss,ss])
p_v_up1=p.up1
p_v_low1=p.low1
smoo_v_up1=s.up1
smoo_v_low1=s.low1
pre=est.p[1:(S-S0-K+1),ss]
smoo=est.s[(1:S0),ss]
loc1=c(1:S)*fit$gap.len
if(K==1){
smooth_pred=c(smoo,pre)
smooth_pred_up=c(smoo_v_up1,p_v_up1)
smooth_pred_low=c(smoo_v_low1,p_v_low1)
}
if(K>1){
betw=rep(NA,K-1)
smooth_pred=c(smoo,betw,pre)
smooth_pred_up=c(smoo_v_up1,betw,p_v_up1)
smooth_pred_low=c(smoo_v_low1,betw,p_v_low1)
}
if(ss %in% (2*(1:rows)-1)){
#lowbound=min(na.omit(smooth_pred_low))#-0.5*diff(range(na.omit(c(smoo_v_low1,p_v_low1))))
#maxbound=max(na.omit(smooth_pred_up))#+0.5*diff(range(na.omit(c(smoo_v_up1,p_v_up1))))
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
lowbound=min(na.omit(smooth_pred))-scal*rang
maxbound=max(na.omit(smooth_pred))+scal*rang
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(beta[.(vary.names[((ss+1)/2)])]~(t)))
}
if(ss %in% (2*rows+1):(2*rows+q2)){
#lowbound=min(na.omit(smooth_pred_low))#-1.5*max(diff(range(na.omit(c(smoo_v_low1,p_v_low1)))),0.2)
#maxbound=max(na.omit(smooth_pred_up))#+1.5*max(diff(range(na.omit(c(smoo_v_up1,p_v_up1)))),0.2)
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
#rang=maxbound-lowbound
lowbound=min(na.omit(smooth_pred))-scal*rang
maxbound=max(na.omit(smooth_pred))+scal*rang
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(alpha[.(con.names[ss-2*rows])]))
}
lines(loc1,smooth_pred,lty=1)
lines(loc1,smooth_pred_up,lty=3)
lines(loc1,smooth_pred_low,lty=3)
vv=fit$gap.len*training_samp
abline(v=vv,col="grey",lty=2)
#legend('bottomright',c("predict","filter","smooth"),lty=c(1,1,1),col=c("blue","black","red"),text.width=0.15,seg.len=0.5)
}
}
if(initial.fit==T){
est.s=fit$smooth
est.var.s=fit$smooth.var
for(ss in plot.index){
s.up1=est.s[,ss]+2*sqrt(est.var.s[,ss,ss])
s.low1=est.s[,ss]-2*sqrt(est.var.s[,ss,ss])
smoo_v_up1=s.up1
smoo_v_low1=s.low1
smoo=est.s[,ss]
loc1=c(1:S0)*fit$gap.len
smooth_pred=c(smoo)
smooth_pred_up=c(smoo_v_up1)
smooth_pred_low=c(smoo_v_low1)
if(ss %in% (2*(1:rows)-1)){
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
#rang=maxbound-lowbound
lowbound=min(na.omit(smooth_pred))-scal*rang
maxbound=max(na.omit(smooth_pred))+scal*rang
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(beta[.(vary.names[((ss+1)/2)])]~(t)))
}
if(ss %in% (2*rows+1):(2*rows+q2)){
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
#rang=maxbound-lowbound
lowbound=min(na.omit(smooth_pred))-scal*rang
maxbound=max(na.omit(smooth_pred))+scal*rang
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(alpha[.(con.names[ss-2*rows])]))
}
lines(loc1,smooth_pred,lty=1)
lines(loc1,smooth_pred_up,lty=3)
lines(loc1,smooth_pred_low,lty=3)
vv=fit$gap.len*training_samp
abline(v=vv,col="grey",lty=2)
#legend('bottomright',c("predict","filter","smooth"),lty=c(1,1,1),col=c("blue","black","red"),text.width=0.15,seg.len=0.5)
}
}
}
if(S0==S){
est.s=fit$smooth
est.var.s=fit$smooth.var
for(ss in plot.index){
s.up1=est.s[,ss]+2*sqrt(est.var.s[,ss,ss])
s.low1=est.s[,ss]-2*sqrt(est.var.s[,ss,ss])
smoo_v_up1=s.up1
smoo_v_low1=s.low1
smoo=est.s[,ss]
loc1=c(1:S)*fit$gap.len
smooth_pred=c(smoo)
smooth_pred_up=c(smoo_v_up1)
smooth_pred_low=c(smoo_v_low1)
if(ss %in% (2*(1:rows)-1)){
#lowbound=min(na.omit(smooth_pred_low))#-0.5*diff(range(na.omit(c(smoo_v_low1,p_v_low1))))
#maxbound=max()#+0.5*diff(range(na.omit(c(smoo_v_up1,p_v_up1))))
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
lowbound=min(smooth_pred)-rang*scal
maxbound=max(smooth_pred)+rang*scal
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(beta[.(vary.names[((ss+1)/2)])]~(t)))
}
if(ss %in% (2*rows+1):(2*rows+q2)){
#lowbound=min(na.omit(smooth_pred_low))#-1.5*max(diff(range(na.omit(c(smoo_v_low1,p_v_low1)))),0.2)
#maxbound=max(na.omit(smooth_pred_up))#+1.5*max(diff(range(na.omit(c(smoo_v_up1,p_v_up1)))),0.2)
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
lowbound=min(na.omit(smooth_pred))-rang*scal
maxbound=max(na.omit(smooth_pred))+rang*scal
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(alpha[.(con.names[ss-2*rows])]))
}
lines(loc1,smooth_pred,lty=1)
lines(loc1,smooth_pred_up,lty=3)
lines(loc1,smooth_pred_low,lty=3)
#vv=fit$gap.len*training_samp
#abline(v=vv,col="grey",lty=2)
#legend('bottomright',c("predict","filter","smooth"),lty=c(1,1,1),col=c("blue","black","red"),text.width=0.15,seg.len=0.5)
}
}
}
jiakunj/DLSSM documentation built on Dec. 31, 2022, 2:18 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions DLSSM.plot DLSSM DLSSM.valid DLSSM.smooth DLSSM.filter DLSSM.predict DLSSM.init Batched
Documented in DLSSM.filter DLSSM.init DLSSM.plot DLSSM.predict DLSSM.smooth DLSSM.valid
#' Combine data into Batched data
#' @author Jiakun Jiang, Wei Yang, Wensheng Guo
#' @description The time domain of observation will first be standardized into [0,1]. Then [0,1] will be divided into S equally spaced intervals as described in Jiang et al.(2021, Biometrics).
#' Then those intervals slice the dataset to S batches of data.
#' @param formula An object of class "formula" (or one that can be coerced to that class): a symbolic description of response and covariates in the model.
#' @param data Dataset matrix containing the observations (one rows is a sample).
#' @param time The time variable in the dataset. The varying coefficient functions are assumed to be smooth functions of this variable.
#' @param S Number of batches
#' @export
#' @return
#' \tabular{ll}{
#' \code{batched} \tab List of batched data, the element of list is matrix with each row representing a sample \cr
#' \tab \cr
#' \code{gap.len} \tab interval length 1/S \cr
#' }
Batched<-function(formula,data,time,S){
if(any(is.na(data))==T) { stop("The dataset have missing data") }
var_names=all.vars(formula)
if(all(c(var_names,time)%in%colnames(data))==F) {stop("There exist undefined variables in formula")}
y=data[,var_names[1]]
x=data[,var_names[-1]]
t=data[,time]
n=dim(x)[1] #sample size
q=dim(x)[2] #dimension of covariates
if(n!=length(t)|n!=length(y)){
stop("The dimension of data does not match.")
}
if(as.integer(abs(S))!=S){
stop("S need to be a positive integer")
}
t=t/max(t)
oder=order(t)
t.order=t[oder]
x.order=x[oder,]
y.order=y[oder]
# create batch
t.batch=list() #Observational time-points
x.batch=list() #Covariates corresponding to varying coefficients
y.batch=list() #Binary outcome
data.all=list()
Data.check=rep(NA,S)
for(i in 1:S){
sel=t.order>(i-1)/S & t.order<=i/S
t.batch[[i]]=t.order[sel]
x.batch[[i]]=x.order[sel,]
y.batch[[i]]=as.numeric(y.order[sel])
Data.check[i]=sum(sel)
mmm=cbind(y.batch[[i]],x.batch[[i]],t.batch[[i]])
colnames(mmm)=c(var_names,time)
data.all[[i]]=mmm
}
if(any(Data.check==0)){
stop("There exist at least one interval with no data")
}
invisible(list(batched=data.all,x.batch=x.batch,y.batch=y.batch,t.batch=t.batch,gap.len=1/S,S=S,time=time,formula=formula))
#invisible(list(x.batch=x.batch,y.batch=y.batch,t.batch=t.batch,gap.len=1/S))
}
#' Initial model fitting
#' @author Jiakun Jiang, Wei Yang and Wensheng Guo
#' @description This function is for tuning smoothing parameters using training data. The likelihood was calculated by Kalman Filter and maximized to estimate the smoothing parameters.
#' For the given smoothing parameters, the model coefficients can
#' be efficiently estimated using a Kalman filtering algorithm.
#' @param data.batched A object generated by function Data.batched()
#' @param S0 Number of batches of data to be used as training dataset
#' @param vary.effects The names of variables in the dataset assumed to have a time-varying regression effect on the outcome.
#' @param autotune T/F indicates whether or not the automatic tuning procedure described in Jiang et al. (2021) should be applied. Default is true.
#' @param Lambda Specify smoothing parameters if autotune=F
#' @import Matrix
#' @export
#' @importFrom "graphics" "abline" "lines" "par"
#' @importFrom "stats" "na.omit" "optim" "optimize"
#' @return
#' \tabular{ll}{
#' \code{Lambda:} \tab smoothing parameters \cr
#' \tab \cr
#' \code{Smooth:} \tab smoothed state vector \cr
#' \tab \cr
#' \code{Smooth.var:} \tab covariance of smoothed state vector in Smooth. \cr
#' }
#' @examples
#' \donttest{
#' set.seed(321)
#' n=8000
#' beta0=function(t) 0.1*t-1
#' beta1=function(t) cos(2*pi*t)
#' beta2=function(t) sin(2*pi*t)
#' alph1=alph2=1
#' x=matrix(runif(n*4,min=-4,max=4),nrow=n,ncol=4)
#' t=sort(runif(n))
#' coef=cbind(beta0(t),beta1(t),beta2(t),rep(alph1,n),rep(alph2,n))
#' covar=cbind(rep(1,n),x)
#' linear=apply(coef*covar,1,sum)
#' prob=exp(linear)/(1+exp(linear))
#' y=as.numeric(runif(n)<prob)
#' sim.data=cbind(y,x,t)
#' colnames(sim.data)=c("y","x1","x2","x3","x4","t")
#' formula = y~x1+x2+x3+x4
#' # Divide the time domain [0,1] into S=100 equally spaced intervals
#' S=100
#' S0=75
#' data.batched=Batched(formula, data=sim.data, time="t", S)
#'
#' # using first 75 batches as training dataset to tune smoothing parameters
#' fit0=DLSSM.init(data.batched, S0, vary.effects=c("x1","x2"))
#' fit0$Lambda
#' DLSSM.plot(fit0)
#' }
DLSSM.init<-function(data.batched,S0,vary.effects,autotune=TRUE,Lambda=NULL){
formula=data.batched$formula
S=data.batched$S
if(S<S0){stop("Number of batches of training dataset exceed S")}
x.names=all.vars(formula)[-1]
vary=match(vary.effects,x.names)
time=data.batched$time
x.b=data.batched$x.batch
y.b=data.batched$y.batch
t.b=data.batched$t.batch
q=dim(x.b[[1]])[2] #dimension of covariates filter
q1=num.vary=length(vary.effects) #number of covariates with varying coefficient
q2=q-q1 #number of covariates with time-invariant coefficient
if(num.vary>0&num.vary<length(x.names)){
t=list() # observational time-points
X=list() # Covariates corresponding to varying coefficients
XX=list() # Covariates corresponding to constant coefficients
y=list() # binary outcome
for(i in 1:S0){
t[[i]]=t.b[[i]]
X[[i]]=x.b[[i]][,vary]
XX[[i]]=x.b[[i]][,setdiff(1:q,vary)]
y[[i]]=as.numeric(y.b[[i]])
}
}
if(num.vary==0){ # no varying-coefficient
t=list()
XX=list()
y=list()
for(i in 1:S0){
t[[i]]=t.b[[i]]
XX[[i]]=x.b[[i]]
y[[i]]=as.numeric(y.b[[i]])
}
}
if(num.vary==length(x.names)){
t=list() # observational time-points
X=list() # Covariates corresponding to varying coefficients
y=list() # binary outcome
for(i in 1:S0){
t[[i]]=t.b[[i]]
X[[i]]=x.b[[i]][,vary]
y[[i]]=as.numeric(y.b[[i]])
}
}
dim=q1+1 # Number of varying coefficients including intercept
dim.con=q2 # Number of constant coefficients
TT=matrix(0,2*dim+dim.con,2*dim+dim.con) # Transformation matrix in state-space model
QQ=matrix(NA,2,2) # Covariance matrix in state-space model
delta.t=data.batched$gap.len # Batch data be generated equally spaced
TTq1=matrix(0,2*dim,2*dim) # Transformation matrix for varying coefficient
TTq2=diag(rep(1,dim.con)) # Transformation matrix for constant coefficient
for(ss in 1:dim){
TTq1[(2*ss-1):(2*ss),(2*ss-1):(2*ss)]=matrix(c(1, 0, delta.t, 1),2, 2)
}
TT=as.matrix(bdiag(TTq1,TTq2))
QQ=matrix(c((delta.t)^3/3,(delta.t)^2/2,(delta.t)^2/2,delta.t),2,2)
a1=rep(0,2*dim+dim.con) # Initial value
diff=log(100) # Diffuse variance of initial prior distribution
n.train=S0 # number of batches S
#########################################
#Likelihood function to tune smoothing parameters Lambda based on training data
#########################################
if(autotune==TRUE){
likehood.predictive=function(xx) {
Lambda=xx
Q=matrix(0,2*dim+dim.con,2*dim+dim.con)
Q1=matrix(0,2*dim,2*dim)
for(ss in 1:dim){
Q1[(2*ss-1):(2*ss),(2*ss-1):(2*ss)]=exp(Lambda[ss])*QQ
}
Q[1:(2*dim),1:(2*dim)]=Q1
P1 <- exp(diff)*diag(2*dim+dim.con)
Sample.likehood=rep(NA,n.train)
prediction=matrix(NA,n.train,dim*2+dim.con)
prediction.var=array(NA,dim=c(n.train,dim*2+dim.con,dim*2+dim.con))
prediction[1,]=a1
prediction.var[1,,]=P1
filter=matrix(NA,n.train,dim*2+dim.con)
filter.var=array(NA,dim=c(n.train,dim*2+dim.con,dim*2+dim.con))
intial=prediction[1,]
ZZ=matrix(0,2*dim+dim.con,length(t[[1]]))
if(dim>1&dim<(length(x.names)+1)){
ZZ[1,]=1
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[1]])[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[1]])
}
if(dim==1){
ZZ[1,]=1
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[1]])
}
if(dim==(length(x.names)+1)){
ZZ[1,]=1
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[1]])[sss,]
}
#ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[1]])
}
for(inter in 1:1000){
summa=rep(0,2*dim+dim.con)
summa.cov=matrix(0,2*dim+dim.con,2*dim+dim.con)
for(j in 1:length(t[[1]])){
y.hat=exp(ZZ[,j]%*%intial)/(1+exp(ZZ[,j]%*%intial))
summa=summa+as.numeric(y[[1]][j]-y.hat)*ZZ[,j]
summa.cov=summa.cov+as.numeric(y.hat*(1-y.hat))*ZZ[,j]%o%ZZ[,j]
}
SCORE=summa-t(solve(prediction.var[1,,])%*%(intial-matrix(prediction[1,],2*dim+dim.con,1)))
COVMATRIX=-solve(prediction.var[1,,])-summa.cov
recur=intial-(0.2*solve(COVMATRIX)%*%t(SCORE))
befor=intial
intial=recur
if(mean(abs(befor-recur))<0.00001){break}
}
filter[1,]=recur
filter.var[1,,]=solve(-COVMATRIX)
lihood=0
for(j in 1:length(t[[1]])){
y.hatt=exp(ZZ[,j]%*%filter[1,])/(1+exp(ZZ[,j]%*%filter[1,]))
Covmatrix=-solve(prediction.var[1,,])-as.numeric(y.hatt*(1-y.hatt))*ZZ[,j]%o%ZZ[,j]
densi=1/sqrt(det(prediction.var[1,,]))*
exp(-0.5*(filter[1,]-prediction[1,])%*%solve(prediction.var[1,,])%*%(filter[1,]-prediction[1,]))
lihood=lihood+log(2*pi*sqrt(det(solve(-Covmatrix)))*(y.hatt^(y[[1]][j]))*((1-y.hatt)^(1-y[[1]][j]))*densi)
}
Sample.likehood[1]=lihood
for(i in 2:n.train) {
prediction[i,]=TT%*%filter[i-1,]
prediction.var[i,,]=TT%*%filter.var[i-1,,]%*%t(TT)+Q
intial=prediction[i,]
ZZ=matrix(0,dim*2+dim.con,length(t[[i]]))
ZZ[1,]=1 #num.vary>0&num.vary<length(x.names)
if(num.vary>0&num.vary<length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
if(num.vary==0){
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
if(num.vary==length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
}
for(inter in 1:1000){
summa=rep(0,dim*2+dim.con)
summa.cov=matrix(0,dim*2+dim.con,dim*2+dim.con)
for(j in 1:length(t[[i]])){
y.hat=exp(ZZ[,j]%*%intial)/(1+exp(ZZ[,j]%*%intial))
summa=summa+as.numeric(y[[i]][j]-y.hat)*ZZ[,j]
summa.cov=summa.cov+as.numeric(y.hat*(1-y.hat))*ZZ[,j]%o%ZZ[,j]
}
SCORE=summa-t(solve(prediction.var[i,,])%*%(intial-matrix(prediction[i,],2*dim+dim.con,1)))
COVMATRIX=-solve(prediction.var[i,,])-summa.cov
recur=intial-(0.2*solve(COVMATRIX)%*%t(SCORE))
befor=intial
intial=recur
if(max(abs(befor-recur))<0.00001){break}
}
filter[i,]=recur
filter.var[i,,]=solve(-COVMATRIX)
lihood=0
for(j in 1:length(t[[i]])){
y.hatt=exp(ZZ[,j]%*%filter[i,])/(1+exp(ZZ[,j]%*%filter[i,]))
Covmatrix=-solve(prediction.var[i,,])-as.numeric(y.hatt*(1-y.hatt))*ZZ[,j]%o%ZZ[,j]
densi=1/sqrt(det(prediction.var[i,,]))*
exp(-0.5*(filter[i,]-prediction[i,])%*%solve(prediction.var[i,,])%*%(filter[i,]-prediction[i,]))
lihood=lihood+log(2*pi*sqrt(det(solve(-Covmatrix)))*(y.hatt^(y[[i]][j]))*((1-y.hatt)^(1-y[[i]][j]))*densi)
}
Sample.likehood[i]=lihood
}
likelihood=-sum((Sample.likehood))
return(likelihood)
}
if(dim>1){
search.opt=optim(rep(2,dim),likehood.predictive) #Optimize likelihood function with initial value (0,0)
Lambda=exp(search.opt$par)
}
if(dim==1){
search.opt=optimize(likehood.predictive,c(-20,20)) #Optimize likelihood function with initial value (0,0)
Lambda=exp(search.opt$minimum)
}
#Smoothing parameter for intercept
}else{
if(length(Lambda)!=num.vary+1){
stop("The number of smoothing parameters should equal to number of varying coefficients, including varying intercept.")
}
Lambda=Lambda
}
P1<-exp(diff)*diag(dim*2+dim.con) #Diffuse prior
Q=matrix(0,dim*2+dim.con,dim*2+dim.con)
Q1=matrix(0,2*dim,2*dim)
for(ss in 1:dim){
Q1[(2*ss-1):(2*ss),(2*ss-1):(2*ss)]=Lambda[ss]*QQ
}
Q[1:(2*dim),1:(2*dim)]=Q1
#########################################
#Model estimation using all data with Kalman filter
#########################################
prediction=matrix(NA,S0,2*dim+dim.con) # Mean prediction
prediction.var=array(NA,dim=c(S0,2*dim+dim.con,2*dim+dim.con)) # Covariance prediction
prediction[1,]=a1 # Initial mean
prediction.var[1,,]=P1 # Initial covariance
filter=matrix(NA,S0,2*dim+dim.con) # Filtering mean
filter.var=array(NA,dim=c(S0,2*dim+dim.con,2*dim+dim.con)) # Filtering covariance
#Kalman filter on first timepoint
intial=prediction[1,]
ZZ=matrix(0,2*dim+dim.con,length(t[[1]])) #define covariate z in our state space model (5)
ZZ[1,]=1
if(num.vary>0&num.vary<length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[1]])[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[1]])
}
if(num.vary==0){
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[1]])
}
if(num.vary==length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[1]])[sss,]
}
}
for(inter in 1:1000){ # Newton-Raphson iterative algorithm to get filtering estimator
summa=rep(0,2*dim+dim.con)
summa.cov=matrix(0,2*dim+dim.con,2*dim+dim.con)
for(j in 1:length(t[[1]])){
y.hat=exp(ZZ[,j]%*%intial)/(1+exp(ZZ[,j]%*%intial))
summa=summa+as.numeric(y[[1]][j]-y.hat)*ZZ[,j]
summa.cov=summa.cov+as.numeric(y.hat*(1-y.hat))*ZZ[,j]%o%ZZ[,j]
}
SCORE=summa-t(solve(prediction.var[1,,])%*%(intial-matrix(prediction[1,],2*dim+dim.con,1)))
COVMATRIX=-solve(prediction.var[1,,])-summa.cov
recur=intial-(0.2*solve(COVMATRIX)%*%t(SCORE))
befor=intial
intial=recur
if(mean(abs(befor-recur))<0.00001){break}
if(inter==1000){stop("failed converge","\n")}
}
filter[1,]=recur # Filtering mean
filter.var[1,,]=solve(-COVMATRIX) # Filtering covariance
#Kalman filter on timepoints from 2 to S
Prob.est=list()
for(i in 2:S0) {
# Prediction step
prediction[i,]=TT%*%filter[i-1,]
prediction.var[i,,]=TT%*%filter.var[i-1,,]%*%t(TT)+Q
# Filtering step
intial=prediction[i,]
ZZ=matrix(0,dim*2+dim.con,length(t[[i]]))
ZZ[1,]=1
if(num.vary>0&num.vary<length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
if(num.vary==0){
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
if(num.vary==length(x.names)){
for (sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
}
if(i>n.train){
y.hat=exp(t(ZZ)%*%intial)/(1+exp(t(ZZ)%*%intial))
Prob.est[[i-n.train]]=y.hat
}
for(inter in 1:1000){ # Newton-Raphson
summa=rep(0,dim*2+dim.con)
summa.cov=matrix(0,dim*2+dim.con,dim*2+dim.con)
for(j in 1:length(t[[i]])){
y.hat=exp(ZZ[,j]%*%intial)/(1+exp(ZZ[,j]%*%intial))
summa=summa+as.numeric(y[[i]][j]-y.hat)*ZZ[,j]
summa.cov=summa.cov+as.numeric(y.hat*(1-y.hat))*ZZ[,j]%o%ZZ[,j]
}
SCORE=summa-t(solve(prediction.var[i,,])%*%(intial-matrix(prediction[i,],2*dim+dim.con,1)))
COVMATRIX=-solve(prediction.var[i,,])-summa.cov
recur=intial-(0.2*solve(COVMATRIX)%*%t(SCORE))
befor=intial
intial=recur
if(mean(abs(befor-recur))<0.00001){break}
if(inter==1000){stop("Newton-Raphson failed to converge","\n")}
#if(inter==1000){break}
}
filter[i,]=recur
filter.var[i,,]=solve(-COVMATRIX)
}
#Smoothed coefficients on training data
state.smooth=matrix(NA,S0,2*dim+dim.con)
F.t.smooth=array(NA,dim=c(S0,2*dim+dim.con,2*dim+dim.con))
state.smooth[S0,]=filter[S0,]
F.t.smooth[S0,,]=filter.var[S0,,]
for(i in (S0-1):1){
A.i=filter.var[i,,]%*%t(TT)%*%solve(prediction.var[i+1,,])
state.smooth[i,]=filter[i,]+A.i%*%(state.smooth[i+1,]-prediction[i+1,])
F.t.smooth[i,,]=filter.var[i,,]-A.i%*%(prediction.var[i+1,,]-F.t.smooth[i+1,,])%*%t(A.i)
}
Est=list(pred=prediction,pred.var=prediction.var,filter=filter
,filter.var=filter.var,smooth=state.smooth,smooth.var=F.t.smooth
,Lambda=Lambda,vary.effects=vary.effects,q=q,q1=q1,q2=q2
,dim=dim,dim.con=dim.con,num.vary=num.vary,TT=TT,Q=Q,time=time
,Prob.est=as.vector(Prob.est),S=S,S0=S0,gap.len=delta.t,formula=formula,initial.fit=TRUE)
return(Est)
}
DLSSM.predict<-function(fit,K,newx){
vary=fit$vary
num.vary=length(vary)
dim=fit$dim
dim.con=fit$dim.con
if(is.vector(fit$filter)==T){
fil.last=fit$filter
filt.var.last=fit$filter.var
}
else{
fil.last=fit$filter[dim(fit$filter)[1],]
filt.var.last=fit$filter.var[dim(fit$filter)[1],,]
}
coef.prediction=matrix(NA,K,2*dim+dim.con) # Mean prediction
coef.prediction.var=array(NA,dim=c(K,2*dim+dim.con,2*dim+dim.con)) # Covariance prediction
for(pred in 1:K){
TT=fit$TT
Q=fit$Q
coef.prediction[pred,]=TT%*%fil.last
coef.prediction.var[pred,,]=TT%*%filt.var.last%*%t(TT)+Q
}
coef.pred.onlyK=coef.prediction[K,] # only k-step
coef.pred.onlyK.var=coef.prediction.var[K,,]
if(is.null(newx)==TRUE){prob.pred=NULL}
if(is.null(newx)==FALSE){
if(is.vector(newx)==T) {newx=matrix(newx,nrow=1)}
prediction=coef.pred.onlyK
x.names=all.vars(fit$formula)[-1]
if(all(c(x.names)%in%colnames(newx))==F) {stop("There exist variable in model but not in data")}
newx=newx[,x.names]
vary.effects=fit$vary.effects
vary=match(vary.effects,x.names[-1])
q=dim(newx)[2]
newx=newx[,c(vary,setdiff(1:q,vary))]
ZZ=matrix(0,dim(newx)[1],dim*2+dim.con)
ZZ[,1]=1
for(sss in 1:(dim-1)){
ZZ[,2*sss+1]=newx[,sss]
}
ZZ[,(2*dim+1):(2*dim+dim.con)]=newx[,dim:fit$q]
prob.pred=exp(ZZ%*%prediction)/(1+exp(ZZ%*%(prediction)))
}
return(list(coef.pred=coef.pred.onlyK,coef.pred.var=coef.pred.onlyK.var,prob.pred=prob.pred))
}
DLSSM.filter<-function(fit,newdata){
formula=fit$formula
x.names=all.vars(formula)[-1]
time=fit$time
if(all(c(x.names,time)%in%colnames(newdata))==F) {stop("There exist variable in model but not in data")}
newx=newdata[,x.names]
newy=newdata[,all.vars(formula)[1]]
vary.effects=fit$vary.effects
vary=match(vary.effects,x.names)
q=dim(newx)[2]
newx=newx[,c(vary,setdiff(1:q,vary))]
num.vary=length(vary.effects)
dim=fit$dim
dim.con=fit$dim.con
TT=fit$TT
Q=fit$Q
if(is.vector(fit$filter)==T){
fil.last=fit$filter
filt.var.last=fit$filter.var
}
if(is.vector(fit$filter)==F){
fil.last=fit$filter[dim(fit$filter)[1],]
filt.var.last=fit$filter.var[dim(fit$filter)[1],,]
}
predict=TT%*%fil.last
predict.var=TT%*%filt.var.last%*%t(TT)+Q
# Filtering step
intial=predict
ZZ=matrix(0,2*dim+dim.con,length(newy))
ZZ[1,]=1#num.vary>0&num.vary<length(x.names)
if(num.vary>0&num.vary<length(x.names)){
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(newx)[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(newx)[dim:q,]
}
if(num.vary==0){
ZZ[(2*dim+1):(2*dim+dim.con),]=t(newx)
}
if(num.vary==length(x.names)){
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(newx)[sss,]
}
#ZZ[(2*dim+1):(2*dim+dim.con),]=t(newx)[dim:q,]
}
for(inter in 1:1000){ # Newton-Raphson
summa=rep(0,dim*2+dim.con)
summa.cov=matrix(0,dim*2+dim.con,dim*2+dim.con)
for(j in 1:length(newy)){
y.hat=exp(ZZ[,j]%*%intial)/(1+exp(ZZ[,j]%*%intial))
summa=summa+as.numeric(newy[j]-y.hat)*ZZ[,j]
summa.cov=summa.cov+as.numeric(y.hat*(1-y.hat))*ZZ[,j]%o%ZZ[,j]
}
SCORE=summa-t(solve(predict.var)%*%(intial-matrix(predict,2*dim+dim.con,1)))
COVMATRIX=-solve(predict.var)-summa.cov
recur=intial-(0.2*solve(COVMATRIX)%*%t(SCORE))
befor=intial
intial=recur
if(mean(abs(befor-recur))<0.00001){break}
if(inter==1000){stop("failed converge","\n")}
}
filter=as.vector(recur)
filter.var=solve(-COVMATRIX)
Lambda=fit$Lambda
q=fit$q
q1=fit$q1
Est=list(filter=filter,filter.var=filter.var,Lambda=Lambda,vary=vary,q=q,q1=q1,dim=dim,dim.con=dim.con,TT=TT,Q=Q,time=time,formula=formula,vary.effects=vary.effects)
return(Est)
}
DLSSM.smooth<-function(fit,prediction,prediction.var,filter,filter.var){
dim=fit$dim
dim.con=fit$dim.con
TT=fit$TT
n.train=dim(filter)[1]
state.smooth=matrix(NA,n.train,2*dim+dim.con)
F.t.smooth=array(NA,dim=c(n.train,2*dim+dim.con,2*dim+dim.con))
state.smooth[n.train,]=filter[n.train,]
F.t.smooth[n.train,,]=filter.var[n.train,,]
for(i in (n.train-1):1){
A.i=filter.var[i,,]%*%t(TT)%*%solve(prediction.var[i+1,,])
state.smooth[i,]=filter[i,]+A.i%*%(state.smooth[i+1,]-prediction[i+1,])
F.t.smooth[i,,]=filter.var[i,,]-A.i%*%(prediction.var[i+1,,]-F.t.smooth[i+1,,])%*%t(A.i)
}
return(list(smooth=state.smooth,smooth.var=F.t.smooth))
}
#' Dynamical prediction on validation dataset
#' @author Jiakun Jiang, Wei Yang and Wensheng Guo
#' @description After we have fitted initial model, we can do validation. It is iteratively doing K-steps ahead prediction and model updating (filtering)
#' when a new batch of data becomes available.
#' The validation include K-steps ahead prediction of state vector and probabilities on validation interval.
#' @param fit0 Initial fitted model
#' @param data.batched Batched dataset generated by function Batched()
#' @param K Number of steps for ahead prediction
#' @import Matrix
#' @export
#' @details The argument fit could be object of DLSSM or DLSSM.init.
#' @importFrom "graphics" "abline" "lines" "par"
#' @importFrom "stats" "na.omit" "optim" "optimize"
#' @return
#' \tabular{ll}{
#' \code{pred.K:} \tab K-steps ahead predicted coefficients \cr
#' \tab \cr
#' \code{pred.var.K:} \tab covariance of K-steps ahead predicted coefficients \cr
#' \tab \cr
#' \code{pred.prob.K:} \tab K-steps ahead predicted probabilities \cr
#' }
#' @examples
#'\donttest{
#' set.seed(321)
#' n=8000
#' beta0=function(t) 0.1*t-1
#' beta1=function(t) cos(2*pi*t)
#' beta2=function(t) sin(2*pi*t)
#' alph1=alph2=1
#' x=matrix(runif(n*4,min=-4,max=4),nrow=n,ncol=4)
#' t=sort(runif(n))
#' coef=cbind(beta0(t),beta1(t),beta2(t),rep(alph1,n),rep(alph2,n))
#' covar=cbind(rep(1,n),x)
#' linear=apply(coef*covar,1,sum)
#' prob=exp(linear)/(1+exp(linear))
#' y=as.numeric(runif(n)<prob)
#' sim.data=cbind(y,x,t)
#' colnames(sim.data)=c("y","x1","x2","x3","x4","t")
#' formula = y~x1+x2+x3+x4
#' # Divide the time domain [0,1] into S=100 equally spaced intervals
#' S=100
#' S0=75
#' data.batched=Batched(formula, data=sim.data, time="t", S)
#'
#' # using first 75 batches as training dataset to tune smoothing parameters
#' fit0=DLSSM.init(data.batched, S0, vary.effects=c("x1","x2"))
#' fit0$Lambda
#'
#' #After initial model fitting on training data, we move to dynamic prediction
#' fit=DLSSM.valid(fit0, data.batched, K=1)
#' DLSSM.plot(fit)
#' }
DLSSM.valid<-function(fit0,data.batched,K){
S=fit0$S
S0=fit0$S0
#if(S0+K>S){stop("S0+K>S happen, the predicted horizon exceed the whole dataset")}
formula=fit0$formula
vary.effects=fit0$vary.effects
TT=fit0$TT
Q=fit0$Q
x.names=all.vars(formula)[-1]
vary=match(vary.effects,x.names)
x.b=data.batched$x.batch
y.b=data.batched$y.batch
t.b=data.batched$t.batch
q=dim(x.b[[1]])[2] #dimension of covariates filter
q1=num.vary=length(vary.effects) #number of covariates with varying coefficient
q2=q-q1 #number of covariates with time-invariant coefficient
if(num.vary>0&num.vary<length(x.names)){
t=list() # observational time-points
X=list() # Covariates corresponding to varying coefficients
XX=list() # Covariates corresponding to constant coefficients
y=list() # binary outcome
for(i in 1:S){
t[[i]]=t.b[[i]]
X[[i]]=x.b[[i]][,vary]
XX[[i]]=x.b[[i]][,setdiff(1:q,vary)]
y[[i]]=as.numeric(y.b[[i]])
}
}
if(num.vary==0){ # no varying-coefficient covariates
t=list()
XX=list()
y=list()
for(i in 1:S){
t[[i]]=t.b[[i]]
XX[[i]]=x.b[[i]]
y[[i]]=as.numeric(y.b[[i]])
}
}
if(num.vary==length(x.names)){
t=list() # observational time-points
X=list() # Covariates corresponding to varying coefficients
#XX=list() # Covariates corresponding to constant coefficients
y=list() # binary outcome
for(i in 1:S){
t[[i]]=t.b[[i]]
X[[i]]=x.b[[i]][,vary]
#XX[[i]]=x.b[[i]][,setdiff(1:q,vary)]
y[[i]]=as.numeric(y.b[[i]])
}
}
dim=fit0$dim
dim.con=fit0$dim.con
dimen=2*dim+dim.con
Prediction=matrix(NA,S,dimen) # 3 varying-coefficients and 2 constants coefficients
Prediction.var=array(NA,dim=c(S,dimen,dimen))
Filtering=matrix(NA,S,dimen)
Filtering.var=array(NA,dim=c(S,dimen,dimen))
Prediction[1:S0,]=fit0$pred
Prediction.var[1:S0,,]=fit0$pred.var
Filtering[1:S0,]=fit0$filter
Filtering.var[1:S0,,]=fit0$filter.var
fit.last=fit0
for(i in (S0+1):S){
pred=DLSSM.predict(fit.last,newx=NULL,K=1)
Prediction[i,]=pred$coef.pred
Prediction.var[i,,]=pred$coef.pred.var
fit.last=DLSSM.filter(fit.last,data.batched$batched[[i]])
Filtering[i,]=fit.last$filter
Filtering.var[i,,]=fit.last$filter.var
}
prediction.K=matrix(NA,S,2*dim+dim.con) #for K steps ahead coef prediction
prediction.var.K=array(NA,dim=c(S,2*dim+dim.con,2*dim+dim.con))
Prob.est.K=list() # K steps ahead prediction
for(i in (S0+K):S){
K.S.P=matrix(NA,K,2*dim+dim.con)
K.S.P.var=array(NA,dim=c(K,2*dim+dim.con,2*dim+dim.con))
K.S.P[1,]=TT%*%Filtering[i-K,]
K.S.P.var[1,,]=TT%*%Filtering.var[i-K,,]%*%t(TT)+Q
if(K>1){
for(kp in 2:K){
K.S.P[kp,]=TT%*%K.S.P[kp-1,]
K.S.P.var[kp,,]=TT%*%K.S.P.var[kp-1,,]%*%t(TT)+Q
}
}
prediction.K[i,]=K.S.P[K,]
prediction.var.K[i,,]=K.S.P.var[K,,]
#K-steps ahead prediction of probability.
y.hat.k=list()
intial=prediction.K[i,]
ZZ=matrix(0,dim*2+dim.con,length(t[[i]]))
if(num.vary==q){
ZZ[1,]=1
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
}
if(num.vary>0&num.vary<q){
ZZ[1,]=1
for(sss in 1:(dim-1)) {
ZZ[2*sss+1,]=t(X[[i]])[sss,]
}
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
if(num.vary==0){
ZZ[1,]=1
ZZ[(2*dim+1):(2*dim+dim.con),]=t(XX[[i]])
}
y.hat.k=exp(t(ZZ)%*%intial)/(1+exp(t(ZZ)%*%intial))
Prob.est.K[[i]]=y.hat.k
}
Smoothed=DLSSM.smooth(fit.last,Prediction,Prediction.var,Filtering,Filtering.var)
Est=list(pred=Prediction,pred.var=Prediction.var,filter=Filtering
,filter.var=Filtering.var,smooth=Smoothed$smooth,smooth.var=Smoothed$smooth.var
,pred.K=prediction.K[(S0+1):S,],pred.var.K=prediction.var.K[(S0+1):S,,],pred.prob.K=as.vector(Prob.est.K),Lambda=fit0$Lambda,vary.effects=vary.effects,q=q,q1=q1,q2=q2
,dim=dim,dim.con=dim.con,TT=TT,Q=Q
,S=S,S0=S0,gap.len=fit0$gap.len,K=K,formula=formula,initial.fit=FALSE,fit0=fit0)
return(Est)
}
#' Combine model training and validation in a integrated function
#' @author Jiakun Jiang, Wei Yang and Wensheng Guo
#' @description This combine model training and validation in a integrated automatic function DLSSM().
#' @param data.batched A object generated by function Data.batched()
#' @param S0 Number of batches of data to be used as training dataset
#' @param vary.effects The names of variables in the dataset assumed to have a time-varying regression effect on the outcome.
#' @param autotune T/F indicates whether or not the automatic tuning procedure desribed in Jiakun et al. (2021) should be applied. Default is true.
#' @param Lambda Specify smoothing parameters if autotune=F
#' @param K Number of steps for ahead prediction
#' @import Matrix
#' @importFrom "graphics" "abline" "lines" "par"
#' @importFrom "stats" "na.omit" "optim" "optimize"
#' @export
#' @return
#' \tabular{ll}{
#' \code{Lambda:} \tab smoothing parameters \cr
#' \tab \cr
#' \code{Smooth:} \tab smoothed state vector \cr
#' \tab \cr
#' \code{Smooth.var:} \tab covariance of smoothed state vector in Smooth. \cr
#' }
#' @examples
#' \donttest{
#' set.seed(321)
#' n=8000
#' beta0=function(t) 0.1*t-1
#' beta1=function(t) cos(2*pi*t)
#' beta2=function(t) sin(2*pi*t)
#' alph1=alph2=1
#' x=matrix(runif(n*4,min=-4,max=4),nrow=n,ncol=4)
#' t=sort(runif(n))
#' coef=cbind(beta0(t),beta1(t),beta2(t),rep(alph1,n),rep(alph2,n))
#' covar=cbind(rep(1,n),x)
#' linear=apply(coef*covar,1,sum)
#' prob=exp(linear)/(1+exp(linear))
#' y=as.numeric(runif(n)<prob)
#' sim.data=cbind(y,x,t)
#' colnames(sim.data)=c("y","x1","x2","x3","x4","t")
#' formula = y~x1+x2+x3+x4
#' # Divide the time domain [0,1] into S=100 equally spaced intervals
#' S=100
#' S0=75
#' data.batched=Batched(formula, data=sim.data, time="t", S)
#'
#'# Take first S0=75 batches as training data, remaining S-S0=25 batches of data as validation data.
#' fit1=DLSSM(data.batched, S0, vary.effects=c("x1","x2"), autotune=TRUE, Lambda=NULL, K=1)
#' DLSSM.plot(fit1)
#' fit2=DLSSM(data.batched, S0, vary.effects=c("x1","x2"), autotune=TRUE, Lambda=NULL, K=2)
#' DLSSM.plot(fit2)
#' }
DLSSM<-function(data.batched, S0, vary.effects, autotune=TRUE, Lambda=NULL, K){
fit0=DLSSM.init(data.batched,S0,vary.effects=vary.effects,autotune,Lambda)
fit=DLSSM.valid(fit0, data.batched, K)
Est=list(pred=fit0$pred,pred.var=fit0$pred.var,filter=fit0$filter
,filter.var=fit0$filter.var,smooth=fit0$smooth,smooth.var=fit0$smooth.var
,pred.K=fit$pred.K,pred.var.K=fit$pred.var.K,pred.prob.K=fit$pred.prob.K,Lambda=fit0$Lambda,vary.effects=vary.effects,q=fit$q,q1=fit$q1,q2=fit$q2
,dim=fit$dim,dim.con=fit$dim.con,TT=fit$TT,Q=fit$Q
,S=fit$S,S0=S0,gap.len=fit0$gap.len,K=K,formula=fit$formula,initial.fit=FALSE,fit0=fit0)
return(Est)
}
#' Plot coefficients
#' @author Jiakun Jiang, Wei Yang and Wensheng Guo
#' @description Plot smoothed coefficients in the training part and predicted coefficients in validation part, the two parts are divided by vertical dash line.
#' @param fit fitted object
#' @importFrom "graphics" "abline" "lines" "par"
#' @importFrom "stats" "na.omit" "optim" "optimize"
#' @export
#' @return Figures
#' @details If argument "fit" is an initial fitted model then only smoothed coefficients part are plotted.
DLSSM.plot<-function(fit){
x.names=all.vars(fit$formula)[-1]
vary.names=c(0,fit$vary.effects)
con.names=setdiff(x.names,vary.names)
S=fit$S
S0=fit$S0
q2=fit$q2
K=fit$K
initial.fit=fit$initial.fit
rows=length(fit$vary)+1
training_samp=fit$S0
cc=rows#+q2
c1=ceiling(sqrt(cc))
f1=floor(sqrt(cc))
if(cc==c1*f1){numb.plot=c(f1,c1)}
if(c1*f1>cc){numb.plot=c(f1,c1)}
if(c1*f1<cc){numb.plot=c(c1,c1)}
#if(rows>1&rows<length(x.names)+1){
plot.index=2*(1:rows)-1
#}
#if(rows==length(x.names)+1){
# plot.index=2*(1:rows)-1
#}
#if(rows==1){
# plot.index=2*(1:rows)-1
#}
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mfrow=numb.plot,mar=c(4, 4, 1, 1),oma=c(1,1,1,1))
scal=3
if(S0<S){
if(initial.fit==F){
est.p=fit$pred.K
est.var.p=fit$pred.var.K
fit0=fit$fit0
est.s=fit0$smooth
est.var.s=fit0$smooth.var
for(ss in plot.index){
p.up1=est.p[1:(S-S0-K+1),ss]+2*sqrt(est.var.p[1:(S-S0-K+1),ss,ss])
p.low1=est.p[1:(S-S0-K+1),ss]-2*sqrt(est.var.p[1:(S-S0-K+1),ss,ss])
s.up1=est.s[(1:S0),ss]+2*sqrt(est.var.s[(1:S0),ss,ss])
s.low1=est.s[(1:S0),ss]-2*sqrt(est.var.s[(1:S0),ss,ss])
p_v_up1=p.up1
p_v_low1=p.low1
smoo_v_up1=s.up1
smoo_v_low1=s.low1
pre=est.p[1:(S-S0-K+1),ss]
smoo=est.s[(1:S0),ss]
loc1=c(1:S)*fit$gap.len
if(K==1){
smooth_pred=c(smoo,pre)
smooth_pred_up=c(smoo_v_up1,p_v_up1)
smooth_pred_low=c(smoo_v_low1,p_v_low1)
}
if(K>1){
betw=rep(NA,K-1)
smooth_pred=c(smoo,betw,pre)
smooth_pred_up=c(smoo_v_up1,betw,p_v_up1)
smooth_pred_low=c(smoo_v_low1,betw,p_v_low1)
}
if(ss %in% (2*(1:rows)-1)){
#lowbound=min(na.omit(smooth_pred_low))#-0.5*diff(range(na.omit(c(smoo_v_low1,p_v_low1))))
#maxbound=max(na.omit(smooth_pred_up))#+0.5*diff(range(na.omit(c(smoo_v_up1,p_v_up1))))
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
lowbound=min(na.omit(smooth_pred))-scal*rang
maxbound=max(na.omit(smooth_pred))+scal*rang
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(beta[.(vary.names[((ss+1)/2)])]~(t)))
}
if(ss %in% (2*rows+1):(2*rows+q2)){
#lowbound=min(na.omit(smooth_pred_low))#-1.5*max(diff(range(na.omit(c(smoo_v_low1,p_v_low1)))),0.2)
#maxbound=max(na.omit(smooth_pred_up))#+1.5*max(diff(range(na.omit(c(smoo_v_up1,p_v_up1)))),0.2)
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
#rang=maxbound-lowbound
lowbound=min(na.omit(smooth_pred))-scal*rang
maxbound=max(na.omit(smooth_pred))+scal*rang
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(alpha[.(con.names[ss-2*rows])]))
}
lines(loc1,smooth_pred,lty=1)
lines(loc1,smooth_pred_up,lty=3)
lines(loc1,smooth_pred_low,lty=3)
vv=fit$gap.len*training_samp
abline(v=vv,col="grey",lty=2)
#legend('bottomright',c("predict","filter","smooth"),lty=c(1,1,1),col=c("blue","black","red"),text.width=0.15,seg.len=0.5)
}
}
if(initial.fit==T){
est.s=fit$smooth
est.var.s=fit$smooth.var
for(ss in plot.index){
s.up1=est.s[,ss]+2*sqrt(est.var.s[,ss,ss])
s.low1=est.s[,ss]-2*sqrt(est.var.s[,ss,ss])
smoo_v_up1=s.up1
smoo_v_low1=s.low1
smoo=est.s[,ss]
loc1=c(1:S0)*fit$gap.len
smooth_pred=c(smoo)
smooth_pred_up=c(smoo_v_up1)
smooth_pred_low=c(smoo_v_low1)
if(ss %in% (2*(1:rows)-1)){
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
#rang=maxbound-lowbound
lowbound=min(na.omit(smooth_pred))-scal*rang
maxbound=max(na.omit(smooth_pred))+scal*rang
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(beta[.(vary.names[((ss+1)/2)])]~(t)))
}
if(ss %in% (2*rows+1):(2*rows+q2)){
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
#rang=maxbound-lowbound
lowbound=min(na.omit(smooth_pred))-scal*rang
maxbound=max(na.omit(smooth_pred))+scal*rang
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(alpha[.(con.names[ss-2*rows])]))
}
lines(loc1,smooth_pred,lty=1)
lines(loc1,smooth_pred_up,lty=3)
lines(loc1,smooth_pred_low,lty=3)
vv=fit$gap.len*training_samp
abline(v=vv,col="grey",lty=2)
#legend('bottomright',c("predict","filter","smooth"),lty=c(1,1,1),col=c("blue","black","red"),text.width=0.15,seg.len=0.5)
}
}
}
if(S0==S){
est.s=fit$smooth
est.var.s=fit$smooth.var
for(ss in plot.index){
s.up1=est.s[,ss]+2*sqrt(est.var.s[,ss,ss])
s.low1=est.s[,ss]-2*sqrt(est.var.s[,ss,ss])
smoo_v_up1=s.up1
smoo_v_low1=s.low1
smoo=est.s[,ss]
loc1=c(1:S)*fit$gap.len
smooth_pred=c(smoo)
smooth_pred_up=c(smoo_v_up1)
smooth_pred_low=c(smoo_v_low1)
if(ss %in% (2*(1:rows)-1)){
#lowbound=min(na.omit(smooth_pred_low))#-0.5*diff(range(na.omit(c(smoo_v_low1,p_v_low1))))
#maxbound=max()#+0.5*diff(range(na.omit(c(smoo_v_up1,p_v_up1))))
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
lowbound=min(smooth_pred)-rang*scal
maxbound=max(smooth_pred)+rang*scal
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(beta[.(vary.names[((ss+1)/2)])]~(t)))
}
if(ss %in% (2*rows+1):(2*rows+q2)){
#lowbound=min(na.omit(smooth_pred_low))#-1.5*max(diff(range(na.omit(c(smoo_v_low1,p_v_low1)))),0.2)
#maxbound=max(na.omit(smooth_pred_up))#+1.5*max(diff(range(na.omit(c(smoo_v_up1,p_v_up1)))),0.2)
rang=max(na.omit(smooth_pred_up)-na.omit(smooth_pred_low))
lowbound=min(na.omit(smooth_pred))-rang*scal
maxbound=max(na.omit(smooth_pred))+rang*scal
plot(NULL,xlim=c(0,1),ylim=c(lowbound,maxbound),xlab="t",ylab=bquote(alpha[.(con.names[ss-2*rows])]))
}
lines(loc1,smooth_pred,lty=1)
lines(loc1,smooth_pred_up,lty=3)
lines(loc1,smooth_pred_low,lty=3)
#vv=fit$gap.len*training_samp
#abline(v=vv,col="grey",lty=2)
#legend('bottomright',c("predict","filter","smooth"),lty=c(1,1,1),col=c("blue","black","red"),text.width=0.15,seg.len=0.5)
}
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.