Hi all,
I'm very confused! I've been using the same code for many weeks without
any
bother for various covariates. I'm now looking at another covaraite and
whenever I run the code you can see below I get an error message: "Error in
rep(0, nrow(data)) : invalid 'times' argument"
This code works:
# remove 'missing' cases from data #
snearma <- function(data)
{
for(i in 1:nrow(data)){
if(is.na(data$all.sp)[i])
data$all.sp[i]<-0
if(is.na(data$all.cp)[i])
data$all.cp[i]<-0
if(is.na(data$all.scgtc)[i])
data$all.scgtc[i]<-0
if(is.na(data$all.tc)[i])
data$all.tc[i] <- 0
if(is.na(data$all.ta)[i])
data$all.ta[i] <- 0
if(is.na(data$all.aa)[i])
data$all.aa[i] <- 0
if(is.na(data$all.m)[i])
data$all.m[i] <- 0
if(is.na(data$all.otc)[i])
data$all.otc[i] <- 0
if(is.na(data$all.o)[i])
data$all.o[i] <- 0
}
dummy <- rep(0,nrow(data))
for(i in 1:nrow(data)){
if(data$all.sp[i]==0 && data$all.cp[i]==0 &&
data$all.scgtc[i]==0 &&
data$all.tc[i]==0 && data$all.ta[i]==0 && data$all.aa[i]==0
&&
data$all.m[i]==0 & data$all.otc[i]==0 && data$all.o[i]==0)
dummy[i] <- i
}
return(data[-dummy,])
}
# create smaller dataset with missing cases removed #
nmarma <- snearma(nearma)
# create short stratification variable #
nmrpa <- randp(nmarma)
# create censoring variable for the covariate #
stypea <- function(data)
{
for(i in 1:nrow(data)){
if(is.na(data$all.sp)[i])
data$all.sp[i]<-0
if(is.na(data$all.cp)[i])
data$all.cp[i]<-0
if(is.na(data$all.scgtc)[i])
data$all.scgtc[i]<-0
if(is.na(data$all.tc)[i])
data$all.tc[i] <- 0
if(is.na(data$all.ta)[i])
data$all.ta[i] <- 0
if(is.na(data$all.aa)[i])
data$all.aa[i] <- 0
if(is.na(data$all.m)[i])
data$all.m[i] <- 0
if(is.na(data$all.otc)[i])
data$all.otc[i] <- 0
if(is.na(data$all.o)[i])
data$all.o[i] <- 0
}
stype <- rep(0,nrow(data))
for(i in 1:nrow(data)){
if(data$all.type[i]=="P" && data$all.sp[i]>=1
&& data$all.scgtc[i] == 0)
stype[i] <- 1
if(data$all.type[i]=="P" && data$all.cp[i]>=1
&& data$all.scgtc[i] == 0)
stype[i] <- 1
if(data$all.type[i]=="P" && data$all.scgtc[i]>=1)
stype[i] <- 2
if(data$all.type[i]=="P" && data$all.sp[i]==0 &&
data$all.cp[i]==0 &&
data$all.scgtc[i]==0 && data$all.otc[i]==0 && data$all.o[i]==0)
stype[i] <- 2
if(data$all.type[i]=="G" && data$all.tc[i]>=1
&& data$all.m[i]==0 &&
data$all.ta[i]==0 & data$all.aa[i]==0)
stype[i] <- 3
if(data$all.type[i]=="G" && data$all.m[i]>=1 &&
data$all.tc[i]==0)
stype[i] <- 3
if(data$all.type[i]=="G" && data$all.ta[i]>=1
&& data$all.tc[i]==0)
stype[i] <- 3
if(data$all.type[i]=="G" && data$all.aa[i]>=1
&& data$all.tc[i]==0)
stype[i] <- 3
if(data$all.type[i]=="G" && data$all.m[i]>=1 &&
data$all.tc[i]>=1)
stype[i] <- 3
if(data$all.type[i]=="G" && data$all.ta[i]>=1
&& data$all.tc[i]>=1)
stype[i] <- 3
if(data$all.type[i]=="G" && data$all.aa[i]>=1
&& data$all.tc[i]>=1)
stype[i] <- 3
if(data$all.otc[i]>=1)
stype[i] <- 4
if(data$all.o[i]>=1)
stype[i] <- 4
}
return(stype)
}
fita <-
survdiff(Surv(rem.Remtime,rem.Rcens)~stypea(nmarma)+strata(nmrpa),data=nmarma)
fita
lrpvalue=1-pchisq(fita$chisq,3)
xx <-
cuminc(nmarma$rem.Remtime/365,nmarma$rem.Rcens,stypea(nmarma),strata=nmrpa)
plot(xx,curvlab=c("Simple/Complex","SC+2gentc or
2gentc","TC or My/Ab or
My/Ab+gentc","Other"),lty=1,color=c(2:5),xlab="Time from
randomisation
(years)",ylab="Probability of 12-month remission",main="Time
to 12-month
remission",wh=c(2.0,0.4))
text(4,0.5,cex=0.85,paste("Log-rank
test=",round(fita$chisq,3),"p-value=",round(lrpvalue,3)))
whereas this doesn't:
par <- function(data)
{
dummy <- rep(0,nrow(data))
for(i in 1:nrow(data)){
if(is.na(data$all.frontlob)[i] && is.na(data$all.templob)[i]
&&
is.na(data$all.parlob)[i]
&& is.na(data$all.occlob)[i] && is.na(data$all.notspec)[i])
dummy[i]<- i
}
for(i in 1:nrow(data)){
if(is.na(data$all.frontlob)[i])
data$all.frontlob[i] <- "N"
if(is.na(data$all.templob)[i])
data$all.templob[i] <- "N"
if(is.na(data$all.parlob)[i])
data$all.parlob[i] <- "N"
if(is.na(data$all.occlob)[i])
data$all.occlob[i] <- "N"
if(is.na(data$all.notspec)[i])
data$all.notspec[i] <- "N"
}
for(i in 1:nrow(data)){
if(data$all.frontlob[i]=="N" &&
data$all.templob[i]=="N" &&
data$all.parlob[i]=="N" && data$all.occlob[i]=="N"
&&
data$all.notspec[i]=="N")
dummy[i] <- i
if(data$all.frontlob[i]=="Y" &&
data$all.parlob[i]=="Y")
dummy[i] <- i
}
return(data[-dummy,])
}
shortpar <- par(nearma)
shortrpa <- randp(shortpar)
lobe <- function(data)
{
for(i in 1:nrow(data)){
if(is.na(data$all.frontlob)[i])
data$all.frontlob[i] <- "N"
if(is.na(data$all.templob)[i])
data$all.templob[i] <- "N"
if(is.na(data$all.occlob)[i])
data$all.occlob[i] <- "N"
if(is.na(data$all.parlob)[i])
data$all.parlob[i] <- "N"
if(is.na(data$all.notspec)[i])
data$all.notspec[i] <- "N"
}
lobe <- rep(0,nrow(data))
for(i in 1:nrow(data)){
if(data$all.frontlob[i]=="Y")
lobe[i] <- 2
if(data$all.templob[i]=="Y")
lobe[i] <- 1
if(data$all.occlob[i]=="Y")
lobe[i] <- 3
if(data$all.parlob[i]=="Y")
lobe[i] <- 3
if(data$all.notspec[i]=="Y")
lobe[i] <- 3
}
return(lobe)
}
fit1 <-
survdiff(Surv(rem.Remtime,rem.Rcens)~lobe(shortpar)+strata(shortrpa),data=shortpar)
fit1
lrpvalue=1-pchisq(fit1$chisq,2)
xx <-
cuminc(shortpar$rem.Remtime/365,shortpar$rem.Rcens,lobe(shortpar),strata=shortrpa)
plot(xx,curvlab=c("Temporal", "Frontal",
"Par./Occ./Other"), lty=1,
color=c(2:4), xlab="Time from randomisation (years)",
ylab="Probability of
12-month remission")
text(4,0.1,cex=0.85,paste("Log-rank test=",
round(fit1$chisq,3),"p-value=",
round(lrpvalue,3)))
which makes me think it is a problem with the data but I can't work out what
that problem is so I'd appreciate any suggestions anyone has to offer!
Thank you,
Laura
[[alternative HTML version deleted]]
