| Title: | Time-Dependent Prognostic Accuracy with Multiply Evaluated Bio Markers or Scores |
|---|---|
| Description: | Time-dependent Receiver Operating Characteristic curves, Area Under the Curve, and Net Reclassification Indexes for repeated measures. It is based on methods in Barbati and Farcomeni (2017) <doi:10.1007/s10260-017-0410-2>. |
| Authors: | Alessio Farcomeni |
| Maintainer: | Alessio Farcomeni <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.0 |
| Built: | 2024-11-18 06:34:46 UTC |
| Source: | CRAN |
Compute area under the ROC curve
auc(ss)auc(ss)
ss |
Matrix with two columns (1-specificities,
sensitivities). It can be simply the output of |
Area under the ROC curve.
A scalar with the AUC.
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods & Applications, in press
# parameters n=100 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker ro=roc(S2,Ti,delta,u,tt,s,vtimes) auc(ro) ## an unrelated marker ro=roc(S1,Ti,delta,u,tt,s,vtimes) auc(ro)# parameters n=100 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker ro=roc(S2,Ti,delta,u,tt,s,vtimes) auc(ro) ## an unrelated marker ro=roc(S1,Ti,delta,u,tt,s,vtimes) auc(ro)
Boostrap the AUC for significance testing and confidence interval calculation
butstrap(X,etime,status,u=NULL,tt,s,vtimes,auc1,B=50,fc=NULL)butstrap(X,etime,status,u=NULL,tt,s,vtimes,auc1,B=50,fc=NULL)
X |
n by S matrix of longitudinal score/biomarker for i-th subject at j-th occasion (NA if unmeasured) |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and
specificity. Defaults to |
tt |
Upper limit (time-horizon) for evaluation of sensitivity and specificity. |
s |
Scalar number of measurements/visits to use for each subject. s<=S |
vtimes |
S vector with visit times |
auc1 |
AUC for the original data set |
B |
Number of bootstrap replicates. Defaults to |
fc |
Events are defined as fc = 1. Defaults to $I(cup X(t_j)>cutoff)$ |
This function can be used to resample the AUC. The resulting p-value is obtained after assumption that the resampled AUC is Gaussian. Non-parametric confidence interval is obtained as the 2.5 and 97.5 confidence interval is simply given by a Gaussian approximation.
A list with the following elements:
p.value |
(Parametric) p-value for H0: AUC=0.5 |
se |
Standard deviation of the AUC replicates |
ci.np |
Non-parametric 95% confidence interval for AUC |
ci.par |
Parametric 95% confidence interval for AUC |
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods \& Applications, in press
# parameters n=100 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## an unimportant marker ro=roc(S1,Ti,delta,u,tt,s,vtimes) but=butstrap(S1,Ti,delta,u,tt,s,vtimes,ro)# parameters n=100 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## an unimportant marker ro=roc(S1,Ti,delta,u,tt,s,vtimes) but=butstrap(S1,Ti,delta,u,tt,s,vtimes,ro)
Boostrap the AUC for significance testing and confidence interval calculation
butstrap.nri(risk1,risk2,etime,status,u,tt,nri1,wh,B=1000)butstrap.nri(risk1,risk2,etime,status,u,tt,nri1,wh,B=1000)
risk1 |
Baseline risk measurements |
risk2 |
Enhanced risk measurements |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and specificity |
tt |
Upper limit (time-horizon) for evaluation of sensitivity and specificity. |
nri1 |
NRI for the original data set |
wh |
Which NRI to boostrap? |
B |
Number of bootstrap replicates. Defaults to |
This function can be used to resample the NRI. The resulting p-value is obtained after assumption that the resampled NRI is Gaussian. Non-parametric confidence interval is obtained as the 2.5 and 97.5 confidence interval is simply given by a Gaussian approximation.
A list with the following elements:
p.value |
(Parametric) p-value for H0: NRI=0 |
se |
Standard deviation of the NRI replicates |
ci.np |
Non-parametric 95% confidence interval for NRI |
ci.par |
Parametric 95% confidence interval for NRI |
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods \& Applications, in press
# parameters n=25 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=1000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] risk1=apply(S1[,1:s],1,sum) risk1=(risk1-min(risk1))/(max(risk1)-min(risk1)) risk2=apply(S2[,1:s],1,sum) risk2=(risk2-min(risk2))/(max(risk2)-min(risk2)) butstrap.nri(risk1,risk2,Ti,delta,u,tt,nri(risk1,risk2,Ti,delta,u,tt)$nri,wh=1,B=500)# parameters n=25 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=1000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] risk1=apply(S1[,1:s],1,sum) risk1=(risk1-min(risk1))/(max(risk1)-min(risk1)) risk2=apply(S2[,1:s],1,sum) risk2=(risk2-min(risk2))/(max(risk2)-min(risk2)) butstrap.nri(risk1,risk2,Ti,delta,u,tt,nri(risk1,risk2,Ti,delta,u,tt)$nri,wh=1,B=500)
Boostrap the AUC for significance testing and confidence interval calculation
butstrap.s(X,etime,status,u=NULL,tt,s,vtimes,auc1,B=50,fc=NULL)butstrap.s(X,etime,status,u=NULL,tt,s,vtimes,auc1,B=50,fc=NULL)
X |
n by S matrix of longitudinal score/biomarker for i-th subject at j-th occasion (NA if unmeasured) |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and
specificity. Defaults to |
tt |
Upper limit (time-horizon) for evaluation of sensitivity and specificity. |
s |
n vector of number of measurements/visits to use for each subject. all(s<=S) |
vtimes |
S vector with visit times |
auc1 |
AUC for the original data set |
B |
Number of bootstrap replicates. Defaults to |
fc |
Events are defined as fc = 1. Defaults to $I(cup X(t_j)>cutoff)$ |
This function can be used to resample the AUC. The resulting p-value is obtained after assumption that the resampled AUC is Gaussian. Non-parametric confidence interval is obtained as the 2.5 and 97.5 confidence interval is simply given by a Gaussian approximation.
A list with the following elements:
p.value |
(Parametric) p-value for H0: AUC=0.5 |
se |
Standard deviation of the AUC replicates |
ci.np |
Non-parametric 95% confidence interval for AUC |
ci.par |
Parametric 95% confidence interval for AUC |
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods \& Applications, in press
# parameters n=100 tt=3 Tmax=10 u=1.5 s=sample(3,n,replace=TRUE) vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## an unimportant marker ro=roc.s(S1,Ti,delta,u,tt,s,vtimes) but=butstrap.s(S1,Ti,delta,u,tt,s,vtimes,ro)# parameters n=100 tt=3 Tmax=10 u=1.5 s=sample(3,n,replace=TRUE) vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## an unimportant marker ro=roc.s(S1,Ti,delta,u,tt,s,vtimes) but=butstrap.s(S1,Ti,delta,u,tt,s,vtimes,ro)
Compute optimal score for AUC
maxauc(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL)maxauc(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL)
X |
p by n by S array of longitudinal scores/biomarkers for i-th subject at j-th occasion (NA if unmeasured) |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and
specificity. Defaults to |
tt |
Upper limit (time-horizon) for evaluation of sensitivity and specificity. |
s |
Scalar number of measurements/visits to use for each subject. s<=S |
vtimes |
S vector with visit times |
fc |
Events are defined as fc = 1. Defaults to $I(cup X(t_j)>cutoff)$ |
This function can be used to find an optimal linear combination of p scores/biomarkers repeatedly measured over time. The resulting score is optimal as it maximizes the AUC among all possible linear combinations. The first biomarker in array X plays a special role, as by default its coefficient is unitary.
A list with the following elements:
beta |
Beta coefficients for the optimal score |
score |
Optimal score |
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods \& Applications, in press
# parameters n=25 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=500 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## X=array(NA,c(2,nrow(S1),ncol(S1))) X[1,,]=round(S2) #fewer different values, quicker computation X[2,,]=S1 sc=maxauc(X,Ti,delta,u,tt,s,vtimes) # beta coefficients sc$beta # final score (X[1,,]+X[2,,]*sc$beta[1]+...+X[p,,]*sc$beta[p-1]) sc$score# parameters n=25 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=500 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## X=array(NA,c(2,nrow(S1),ncol(S1))) X[1,,]=round(S2) #fewer different values, quicker computation X[2,,]=S1 sc=maxauc(X,Ti,delta,u,tt,s,vtimes) # beta coefficients sc$beta # final score (X[1,,]+X[2,,]*sc$beta[1]+...+X[p,,]*sc$beta[p-1]) sc$score
Compute optimal score for AUC
maxauc.s(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL)maxauc.s(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL)
X |
n by S matrix of longitudinal score/biomarker for i-th subject at j-th occasion (NA if unmeasured) |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and
specificity. Defaults to |
tt |
Upper limit (time-horizon) for evaluation of sensitivity and specificity. |
s |
n vector of number of measurements/visits to use for each subject. all(s<=S) |
vtimes |
S vector with visit times |
fc |
Events are defined as fc = 1. Defaults to $I(cup X(t_j)>cutoff)$ |
This function can be used to find an optimal linear combination of p scores/biomarkers repeatedly measured over time. The resulting score is optimal as it maximizes the AUC among all possible linear combinations. The first biomarker in array X plays a special role, as by default its coefficient is unitary.
A list with the following elements:
beta |
Beta coefficients for the optimal score |
score |
Optimal score |
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods \& Applications, in press
# parameters n=20 tt=3 Tmax=10 u=1.5 s=sample(3,n,replace=TRUE) vtimes=c(0,1,2,5) # generate data ngrid=500 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## X=array(NA,c(2,nrow(S1),ncol(S1))) X[1,,]=round(S2) #fewer different values, quicker computation X[2,,]=S1 sc=maxauc.s(X,Ti,delta,u,tt,s,vtimes) # beta coefficients sc$beta # final score (X[1,,]+X[2,,]*sc$beta[1]+...+X[p,,]*sc$beta[p-1]) sc$score# parameters n=20 tt=3 Tmax=10 u=1.5 s=sample(3,n,replace=TRUE) vtimes=c(0,1,2,5) # generate data ngrid=500 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## X=array(NA,c(2,nrow(S1),ncol(S1))) X[1,,]=round(S2) #fewer different values, quicker computation X[2,,]=S1 sc=maxauc.s(X,Ti,delta,u,tt,s,vtimes) # beta coefficients sc$beta # final score (X[1,,]+X[2,,]*sc$beta[1]+...+X[p,,]*sc$beta[p-1]) sc$score
Compute NRI
nri(risk1, risk2, etime,status,u,tt)nri(risk1, risk2, etime,status,u,tt)
risk1 |
Baseline risk measures |
risk2 |
Enhanced risk measures |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and specificity. |
tt |
Upper limit (time-horizon) for evaluation of sensitivity and specificity. |
This function gives the continuous NRI to compare two risk measures.
A list with the following elements:
nri |
1/2 NRI |
nri.events |
NRI for events |
nri.nonevents |
NRI for non-events |
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods \& Applications, in press
# parameters n=100 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] risk1=apply(S1[,1:s],1,sum) risk1=(risk1-min(risk1))/(max(risk1)-min(risk1)) risk2=apply(S2[,1:s],1,sum) risk2=(risk2-min(risk2))/(max(risk2)-min(risk2)) nri(risk1,risk2,Ti,delta,u,tt)# parameters n=100 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] risk1=apply(S1[,1:s],1,sum) risk1=(risk1-min(risk1))/(max(risk1)-min(risk1)) risk2=apply(S2[,1:s],1,sum) risk2=(risk2-min(risk2))/(max(risk2)-min(risk2)) nri(risk1,risk2,Ti,delta,u,tt)
Compute area under the ROC curve for several values of time horizon
plotAUC(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL,plot=TRUE)plotAUC(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL,plot=TRUE)
X |
n by S matrix of longitudinal score/biomarker for i-th subject at j-th occasion (NA if unmeasured) |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and
specificity. Defaults to |
tt |
A vector of upper limits (time-horizons) for evaluation of sensitivity and specificity. |
s |
Scalar number of measurements/visits to use for each subject. s<=S |
vtimes |
S vector with visit times |
fc |
Events are defined as fc = 1. Defaults to $I(cup X(t_j)>cutoff)$ |
plot |
Do we plot the AUCs? Defaults to |
Area under the ROC curve is computed for each value of the vector
tt. The resulting vector is returned. If plot=TRUE (which is
the default) also a plot of tt vs AUC is displayed.
A vector with AUCs
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods & Applications, in press
# parameters n=25 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=1000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker aucs=plotAUC(S2,Ti,delta,u,seq(2,5,length=5),s,vtimes)# parameters n=25 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=1000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker aucs=plotAUC(S2,Ti,delta,u,seq(2,5,length=5),s,vtimes)
Compute area under the ROC curve for several values of the time horizon
plotAUC.s(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL,plot=TRUE)plotAUC.s(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL,plot=TRUE)
X |
n by S matrix of longitudinal score/biomarker for i-th subject at j-th occasion (NA if unmeasured) |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and
specificity. Defaults to |
tt |
A vector of upper limits (time-horizons) for evaluation of sensitivity and specificity. |
s |
n vector of measurements/visits to use for each subject. all(s<=S) |
vtimes |
S vector with visit times |
fc |
Events are defined as fc = 1. Defaults to $I(cup X(t_j)>cutoff)$ |
plot |
Do we plot the AUCs? Defaults to |
Area under the ROC curve is computed for each value of the vector
tt. The resulting vector is returned. If plot=TRUE (which is
the default) also a plot of tt vs AUC is displayed.
A vector with AUCs
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods & Applications, in press
# parameters n=25 tt=3 Tmax=10 u=1.5 s=sample(3,n,replace=TRUE) vtimes=c(0,1,2,5) # generate data ngrid=1000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker aucs=plotAUC.s(S2,Ti,delta,u,seq(2,5,length=5),s,vtimes)# parameters n=25 tt=3 Tmax=10 u=1.5 s=sample(3,n,replace=TRUE) vtimes=c(0,1,2,5) # generate data ngrid=1000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker aucs=plotAUC.s(S2,Ti,delta,u,seq(2,5,length=5),s,vtimes)
Plot the ROC curve
plotROC(ro, add=FALSE, col=NULL)plotROC(ro, add=FALSE, col=NULL)
ro |
Matrix with two columns (1-specificities,
sensitivities). It can be simply the output of |
add |
If |
col |
Colour for the ROC curve (defaults to red) |
Plots the area under the ROC curve.
A plot or a new line in an open plot.
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods & Applications, in press
# parameters n=100 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker ro=roc(S2,Ti,delta,u,tt,s,vtimes) plotROC(ro) ## an unrelated marker ro=roc(S1,Ti,delta,u,tt,s,vtimes) plotROC(ro)# parameters n=100 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker ro=roc(S2,Ti,delta,u,tt,s,vtimes) plotROC(ro) ## an unrelated marker ro=roc(S1,Ti,delta,u,tt,s,vtimes) plotROC(ro)
Compute ROC curve
roc(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL)roc(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL)
X |
n by S matrix of longitudinal score/biomarker for i-th subject at j-th occasion (NA if unmeasured) |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and
specificity. Defaults to |
tt |
Upper limit (time-horizon) for evaluation of sensitivity and specificity. |
s |
Scalar number of measurements/visits to use for each subject. s<=S |
vtimes |
S vector with visit times |
fc |
Events are defined as fc = 1. Defaults to $I(cup X(t_j)>cutoff)$ |
ROC curve is defined as the curve given by (1-specificities, sensitivities). Here these are obtained for a time-dependent multiply-measured marker are defined as
Se(t,c,s,u) = Pr(f_c(X(t_1),X(t_2),...,X(t_s_i))| u <= T <= t),
and
Sp(t,c,s,u) = 1-Pr(f_c(X(t_1),X(t_2),...,X(t_s_i)) | T > t)
for some fixed f_c, where c is a cutoff. The default for f_c is that a positive diagnosis is given as soon as any measurement among the s considered is above the threshold.
A matrix with the following columns:
1-spec |
1-Specificities |
sens |
Sensitivities |
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods \& Applications, in press
# parameters n=100 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker ro=roc(S2,Ti,delta,u,tt,s,vtimes) plot(ro,type="l",col="red") abline(a=0,b=1) ## an unrelated marker ro=roc(S1,Ti,delta,u,tt,s,vtimes) plot(ro,type="l",col="red") abline(a=0,b=1)# parameters n=100 tt=3 Tmax=10 u=1.5 s=2 vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker ro=roc(S2,Ti,delta,u,tt,s,vtimes) plot(ro,type="l",col="red") abline(a=0,b=1) ## an unrelated marker ro=roc(S1,Ti,delta,u,tt,s,vtimes) plot(ro,type="l",col="red") abline(a=0,b=1)
Compute ROC curve
roc.s(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL)roc.s(X,etime,status,u=NULL,tt,s,vtimes,fc=NULL)
X |
n by S matrix of longitudinal score/biomarker for i-th subject at j-th occasion (NA if unmeasured) |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and
specificity. Defaults to |
tt |
Upper limit (time-horizon) for evaluation of sensitivity and specificity. |
s |
n vector of measurements/visits to use for each subject. all(s<=S) |
vtimes |
S vector with visit times |
fc |
Events are defined as fc = 1. Defaults to $I(cup X(t_j)>cutoff)$ |
ROC curve is defined as the curve given by (1-specificities, sensitivities). Here these are obtained for a time-dependent multiply-measured marker are defined as
Se(t,c,s,u) = Pr(f_c(X(t_1),X(t_2),...,X(t_s_i))| u <= T <= t),
and
Sp(t,c,s,u) = 1-Pr(f_c(X(t_1),X(t_2),...,X(t_s_i)) | T > t)
for some fixed f_c, where c is a cutoff. The default for f_c is that a positive diagnosis is given as soon as any measurement among the s considered is above the threshold.
A matrix with the following columns:
1-spec |
1-Specificities |
sens |
Sensitivities |
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods \& Applications, in press
# parameters n=100 tt=3 Tmax=10 u=1.5 s=sample(3,n,replace=TRUE) vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker ro=roc.s(S2,Ti,delta,u,tt,s,vtimes) plot(ro,type="l",col="red") abline(a=0,b=1) ## an unrelated marker ro=roc.s(S1,Ti,delta,u,tt,s,vtimes) plot(ro,type="l",col="red") abline(a=0,b=1)# parameters n=100 tt=3 Tmax=10 u=1.5 s=sample(3,n,replace=TRUE) vtimes=c(0,1,2,5) # generate data ngrid=5000 ts=seq(0,Tmax,length=ngrid) X2=matrix(rnorm(n*ngrid,0,0.1),n,ngrid) for(i in 1:n) { sa=sample(ngrid/6,1) vals=sample(3,1)-1 X2[i,1:sa[1]]=vals[1]+X2[i,1:sa[1]] X2[i,(sa[1]+1):ngrid]=vals[1]+sample(c(-2,2),1)+X2[i,(sa[1]+1):ngrid] } S1=matrix(sample(4,n,replace=TRUE),n,length(vtimes)) S2=matrix(NA,n,length(vtimes)) S2[,1]=X2[,1] for(j in 2:length(vtimes)) { tm=which.min(abs(ts-vtimes[j])) S2[,j]=X2[,tm]} cens=runif(n) ripart=1-exp(-0.01*apply(exp(X2),1,cumsum)*ts/1:ngrid) Ti=rep(NA,n) for(i in 1:n) { Ti[i]=ts[which.min(abs(ripart[,i]-cens[i]))] } cens=runif(n,0,Tmax*2) delta=ifelse(cens>Ti,1,0) Ti[cens<Ti]=cens[cens<Ti] ## ## an important marker ro=roc.s(S2,Ti,delta,u,tt,s,vtimes) plot(ro,type="l",col="red") abline(a=0,b=1) ## an unrelated marker ro=roc.s(S1,Ti,delta,u,tt,s,vtimes) plot(ro,type="l",col="red") abline(a=0,b=1)
Compute sensitivity and specificity
sensspec(X,etime,status,u=NULL,tt,s,vtimes,cutoff=0,fc=NULL)sensspec(X,etime,status,u=NULL,tt,s,vtimes,cutoff=0,fc=NULL)
X |
n by S matrix of longitudinal score/biomarker for i-th subject at j-th occasion (NA if unmeasured) |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and
specificity. Defaults to |
tt |
Upper limit (time-horizon) for evaluation of sensitivity and specificity. |
s |
Scalar number of measurements/visits to use for each subject. s<=S |
vtimes |
S vector with visit times |
cutoff |
cutoff for definining events. Defaults to |
fc |
Events are defined as fc = 1. Defaults to $I(cup X(t_j)>cutoff)$ |
Sensitivity and specificities for a time-dependent multiply-measured marker are defined as
Se(t,c,s,u) = Pr(f_c(X(t_1),X(t_2),...,X(t_s_i))| u <= T <= t),
and
Sp(t,c,s,u) = 1-Pr(f_c(X(t_1),X(t_2),...,X(t_s_i)) | T > t)
for some fixed f_c, where c is a cutoff. The default for f_c is that a positive diagnosis is given as soon as any measurement among the s considered is above the threshold.
A vector with the following elements:
sens |
Sensitivity at the cutoff |
spec |
Specificity at the cutoff |
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods \& Applications, in press
Compute sensitivity and specificity
sensspec.s(X,etime,status,u=NULL,tt,s,vtimes,cutoff=0,fc=NULL)sensspec.s(X,etime,status,u=NULL,tt,s,vtimes,cutoff=0,fc=NULL)
X |
n by S matrix of longitudinal score/biomarker for i-th subject at j-th occasion (NA if unmeasured) |
etime |
n vector with follow-up times |
status |
n vector with event indicators |
u |
Lower limit for evaluation of sensitivity and
specificity. Defaults to |
tt |
Upper limit (time-horizon) for evaluation of sensitivity and specificity. |
s |
n vector of measurements/visits to use for each subject. all(s<=S) |
vtimes |
S vector with visit times |
cutoff |
cutoff for definining events. Defaults to |
fc |
Events are defined as fc = 1. Defaults to $I(cup X(t_j)>cutoff)$ |
Sensitivity and specificities for a time-dependent multiply-measured marker are defined as
Se(t,c,s,u) = Pr(f_c(X(t_1),X(t_2),...,X(t_s_i))| u <= T <= t),
and
Sp(t,c,s,u) = 1-Pr(f_c(X(t_1),X(t_2),...,X(t_s_i)) | T > t)
for some fixed f_c, where c is a cutoff. The default for f_c is that a positive diagnosis is given as soon as any measurement among the s considered is above the threshold.
A vector with the following elements:
sens |
Sensitivity at the cutoff |
spec |
Specificity at the cutoff |
Alessio Farcomeni [email protected]
Barbati, G. and Farcomeni, A. (2017) Prognostic assessment of repeatedly measured time-dependent biomarkers, with application to dilated cardiomuopathy, Statistical Methods \& Applications, in press