Hi everyone,
I am attempting to copy code taken from the article "Decision curve analysis: a technical note" (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6123195/)
When I try to make the simple model taken from the article, I get this error:
"Error in ntbft(data = dt, outcome = outcome, from = model1, xstart = xstart, :
unused argument (from = model1)"
I have no explanation for why this is happening, I am 99% sure I am copying the code correctly. I would appreciate the help!
This code requires nothing to download from except one package, so if anyone would like to copy the code themselves to see if they get this error, here it is:
set.seed(123)
n<-500
rr<-round(abs(rnorm(n,30,10)))
hr<-round(abs(rnorm(n,90,20)))
crp<-round(abs(rnorm(n,150,80)))
install.packages("dummies")
beta0=-7;betarr=0.05
betahr=0.02;betacrp=0.02
linpred<-cbind(1,rr,hr,crp)%*%
c(beta0,betarr,betahr,betacrp)
pi<-exp(linpred)/(1+exp(linpred))
sepsis.tag<-rbinom(n=n,size=1,prob=pi)
dt<-data.frame(rr,hr,crp,sepsis.tag)
ntbft<-function(data,outcome, frm = NULL,
exterdt=NULL, pred=NULL, xstart=0.01,
xstop=0.99, step=0.01, type="treated"){
pt<-seq(from=xstart,to=xstop,by=step)
lpt<-length(pt)
if(type=="treated") coef<-cbind(rep(1,lpt),rep(0,lpt))
if(type=="untreated") coef<-cbind(rep(0,lpt),rep(1,lpt))
if(type=="overall") coef<-cbind(rep(1,lpt),rep(1,lpt))
if(type=="adapt") coef<-cbind(1-pt,pt)
reponse<-as.vector(t(data[outcome]))
if(is.data.frame(exterdt))response<-as.vector(t(exterdt[outcome]))
event.rate<-mean(response)
nball<-event.rate-(1-event.rate)*pt/(1-pt)
nbnone<-1-event.rate-event.rate*(1-pt)/pt
if(is.null(pred)){
model<-glm(frm,data=data,family=binomial("logit"))
pred<-model$fitted.values
if(is.data.frame(exterdt))
pred<-predict(model,newdata=exterdt,type="response")
}
N<-length(pred)
nbt<-rep(NA,lpt)
nbu<-rep(NA,lpt)
for(t in 1:lpt){
tp<-sum(pred>=pt[t]&response==1)
fp<-sum(pred>=pt[t]&response==0)
fn<-sum(pred<pt[t]&response==1)
tn<-sum(pred<pt[t]&response==0)
nbt[t]<-tp/N-fp/N*(pt[t]/(1-pt[t]))
nbu[t]<-tn/N-fn/N*((1-pt[t])/pt[t])
}
nb<-data.frame(pt)
names(nb)<-"threshold"
nb["all"]<-coef[,1]*nball
nb["none"]<-coef[,2]*nbnone
nb["pred"]<-coef[,1]*nbt+coef[,2]*nbu
return(nb)
}
outcome<-"sepsis.tag"
model1<-sepsis.tag~rr+hr
model2<-sepsis.tag~rr+hr+crp
xstart<-0.01;xstop<-0.99;step<-0.01
type<-"treated"
nb.simple<-ntbft(data=dt, outcome=outcome,
from=model1,xstart=xstart,xstop=xstop,
step=step, type=type)