Nothing
Cumrates <-
function (irate,Bdata)
{ # Calculate MSLT from original data (Bdata) and radix
# irate = 1 estimate Nelson-Aalenniter5
# irate = 2 estimate occurrence-exposure rates
# irate = 3 both
# NOTE: if there is an absorbing state and st.absorb=NULL, then
# calling mvna produces and error: "there are undefined transitions".
# In that case, remove the transition from aborbing state to censoring.
if (missing(irate)) irate=3
z<- check.par (Bdata)
namstates <- attr(Bdata,"param")$namstates
numstates <- length (namstates)
namage <- attr(Bdata,"param")$namage
iagehigh <- attr(Bdata,"param")$iagehigh
# require (mvna)
print (". . . . . . Removing intrastate transitions . . . . . ")
locpat <- locpath(Bdata)
removed <- Remove.intrastate(Bdata)
Bdata <- removed
z<- StateSpace (Bdata)
st.absorb <- z$absorbstates
# ======== Step 1: ESTIMATE TRANSITION RATES ========
print (". . . . . . Estimating rates . . . . .")
Lambda <- 0
Lambda.oe <- 0
M.oe <- 0
namstates2 <- vector (mode="numeric",length=length(namstates))
for (i in 1:length(namstates))
{ namstates2[i] <- grep(namstates[i],namstates)}
namstates2 <- as.character(1:length(namstates))
niter5 <- ifelse (irate==2,1,3)
if (irate %in% c(1,3)) # === mvna estimates transition rates ====
{ print (" ")
print ("Cumrates: calls function Biograph.mvna . . . ")
# convert from months to years
Bdata <- date_b (Bdata=Bdata,format.out="age",covs=NULL)
#Bdata <- CMC.ages(Bdata) # NO COVARIATES INCLUDED
print (Bdata[1:5,])
Dmvna <- Biograph.mvna (Bdata)
print ("Cumrates 88")
# see GLHS_mvna.r
D2 <- Dmvna$D
tra <- Dmvna$par$trans_possible
cens <- Dmvna$cens
#D3 <- data.frame(cbind(D2$id,D2$from,as.character(D2$to),D2$entry,D2$exit,D2$time))
D3 <- data.frame(id=D2$id,from=D2$from,
to=D2$to,
entry=round(D2$entry,2),
exit=round(D2$exit,2))
D3$id <- as.numeric(D3$id)
D3$entry <- as.numeric(D3$entry)
D3$exit <- as.numeric(D3$exit)
D3$exit <- ifelse (D3$exit<=D3$entry,D3$entry+1,D3$exit) # CORRECT if CMC exit = CMC entry
# ------------------------------------------------
print (". . . . . . . . . . . . ")
print (". . . . Running mvna . . . . . . ")
# Absorbing states determined in StateSpace (letter)
if (is.null(st.absorb)) zz2 <- D3 else
{ st.absorb2 <- which (namstates==st.absorb)
zz2 <- subset (D3,D3$from!=st.absorb) }
# NOTE: st.absorb used because name of state is given (nog mumber)
na <- mvna(data=zz2,state.names=namstates,tra=attr(Dmvna$D,"param")$trans_possible,cens.name=Dmvna$cens)
cumh <- predict (na,times=seq(0,iagehigh,by=1))
# print ("cumrates test78")
# see MSLT.mvna.r
# Lambda = cumulative hazard by age, destination, origin
nage = nrow(cumh[[1]])
Lambda <- array (NA,dim=c(nage,numstates,numstates,niter5))
if (length(namage)!=nage)
{ print ("WARNING")
print ("Names of ages groups (namage) not consistent with number of age groups")
print ("Names are inferred. Please check the result. ")
namage <- 0:(nage-1)}
dimnames(Lambda) <- list(age=namage[1:nage],destination=namstates,origin=namstates,variant=c("Expected","Upper","Lower"))
for (i in 1:attr(removed,"param")$ntrans)
{ des <- as.numeric(as.character(attr(removed,"param")$transitions$DES[i]))
or <- as.numeric(as.character(attr(removed,"param")$transitions$OR[i]))
Lambda[,des,or,1] <- cumh[[i]]$na
Lambda[,des,or,2] <- cumh[[i]]$upper
Lambda[,des,or,3] <- cumh[[i]]$lower
}
for (iter in 1:3)
{ for (ix in 1:nage)
{ diag(Lambda[ix,,,iter]) <- - apply(Lambda[ix,,,iter],2,function(x) sum(x,na.rm=TRUE)) # column sum
for (i in 1:numstates) {for (j in 1:numstates)
Lambda[ix,j,i,iter] <- ifelse (is.na(Lambda[ix,j,i,iter]),0,Lambda[ix,j,i,iter]) }
}
}
# Compute age-specific rates
astr <- array (NA,dim=c(nage,numstates,numstates,niter5))
dimnames (astr) <- dimnames (Lambda)
astr[1,,,] <- Lambda[1,,,]
for (ix in 2:nage)
{ astr[ix-1,,,] <- Lambda[ix,,,]-Lambda[(ix-1),,,]
}
astr[nage,,,] <- 0
}
if (irate %in% c(2,3))
{ # === occurrence-exposure rates ====
print ("Computing occurrence-exposure rates",quote=FALSE)
print ("Running Sequences.ind",quote=FALSE)
ist <- Sequences.ind (Bdata$path,namstates)
print ("Running Occup",quote=FALSE)
occup <- Occup (Bdata)
print ("Running Trans",quote=FALSE)
trans <- Trans (Bdata)
print ("Running RateTable",quote=FALSE)
ratetable <- RateTable (Bdata,occup,trans)
pc_ac <- 2 # age-cohort rates ( 1 - age-cohort rates)
print ("Running Rates",quote=FALSE)
M <- Rates.ac(ratetable$Stable)
# NOTE: cum oe rates: Lambda.oe[14,,] = sum of rates to and including age 12
Lambda.oe<- M$Mcum
M.oe <- M$M
# compare Lambda[51,,] and Lambda.oe[50,,]
}
if (!exists("astr")) astr <- NA
if (!exists("Lambda")) Lambda <- NA
if (!exists("cumh")) cumh <- NA
if (!exists("Lambda.oe")) Lambda.oe <- NA
if (!exists("M.oe")) M.oe <- NA
cum <- list(D=removed,
irate=irate,
NeAa = Lambda,
predicted=cumh,
astr = astr,
oeCum = Lambda.oe,
oe=M.oe)
class(cum) <- 'cumrates'
return (cum)
# attr(cum$D$D,"param")$namstates
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.