addregmc <- function(addregmce.list,covardata,fromtime=NULL,estcovar=T) {
  # S-PLUS & R functions for analyses of markov surrviel analyses data with an additive hazard model.
  # Written by:	Harald Fekjaer <harald.fekjar@basalmed.uio.no>, Odd O. Aalen and Ørnulf Borgan
  # See also:	http://www.med.uio.no/imb/stat/addreg/
  # License:	GNU version 2
  # Version:	1.0 Beta 1, © 2001
  #-----------Beginning of routine: -----------------
  
  # We check input data:
  if (!is.list(addregmce.list))
    stop("\"addregmce.list\" variable must be an list of \"addregmce\"-objects.")
  if(!is.data.frame(covardata))
    stop("\"covardata\" should be an \"data.frame\".")
  if(dim(covardata)[1]>1)   
    stop("\"covardata\" shall only include one observation.")
  if (sum(names(covardata)=="Constant")>0) 
    stop("\"covardata\" shall not include \"Constant\".")
  else 
    covardata$Constant <- 1
    
  # We read dimentions of data and perform more chechs:
  nrmce <- length(addregmce.list)
  objinfo <- matrix(NA,nrmce,4)
  dimnames(objinfo) <- list(1:nrmce,c("from","to","nrest","mtime"))
  for (i in 1:nrmce) {
    if (class(addregmce.list[[i]])!="addregmce") {
      errormes <- paste("\"addregmce.list\" element nr",i,"is not an \"addregmce\"-object.")
      stop(errormes)
    }
    objinfo[i,1] <- addregmce.list[[i]]$from
    objinfo[i,2] <- addregmce.list[[i]]$to
    objinfo[i,3] <- sum(addregmce.list[[i]]$addreg.object$nevent>0)
    if (addregmce.list[[i]]$addreg.object$stoped)
      objinfo[i,4] <- addregmce.list[[i]]$addreg.object$stopedtime
    else {
      addregtimes  <- addregmce.list[[i]]$addreg.object$times
      objinfo[i,4] <- addregtimes[length(addregtimes)]
    }
  }
  maxnrstate <- max(objinfo[,1:2])
  nrtrans    <- sum(objinfo[,3])

  # We build "vtrans"-matrix with valied transitons.
  vtrans <- matrix(F,maxnrstate,maxnrstate)
  dimnames(vtrans) <- list(1:maxnrstate,1:maxnrstate)
  diag(vtrans) <- T
  for (i in 1:nrmce) {
    from <- objinfo[i,1]
    to   <- objinfo[i,2]
    if (vtrans[from,to]) {
      errormes <- paste("Transition given in element nr",i,"is allready given in an earlier object.")
      stop(errormes)
    }
    vtrans[from,to] <- T
  }

  # We read "addreg" data and chech imput:  
  increhaz <- matrix(NA,nrtrans,5)
  dimnames(increhaz) <- list(1:nrtrans,c("time","dincre","dvar","from","to"))
  startat <- 1 
  for (mcenr in 1:nrmce) {
    nrincre <- objinfo[mcenr,3]
    newincrenr <- startat:(startat+nrincre-1)
    includeincre <- addregmce.list[[mcenr]]$addreg.object$nevent>0
    increhaz[newincrenr,1] <- addregmce.list[[mcenr]]$addreg.object$times[includeincre]
    increhaz[newincrenr,4] <- rep(addregmce.list[[mcenr]]$from,times=nrincre)
    increhaz[newincrenr,5] <- rep(addregmce.list[[mcenr]]$to,times=nrincre)
    startat <- startat+nrincre
    
    # We match covarnames:
    addregincre   <- addregmce.list[[mcenr]]$addreg.object$increments
    if (estcovar)
      addregincrev  <- addregmce.list[[mcenr]]$addreg.object$increcovar
    namesaddreg  <- dimnames(addregincre)[[2]]
    nraddregames <- length(namesaddreg)
    ncovardata    <- names(covardata)
    nrcovardata   <- length(covardata)
    mcovar        <- rep(NA,nraddregames)
    for (mnr in 1:nraddregames) {
      if (sum(namesaddreg[mnr]==namesaddreg)>1) {
        errorm <- paste("Dublicated covarnames in \"addreg.mce\"-object nr",mcenr)
        stop(errorm)
      }
      matches <- namesaddreg[mnr]==ncovardata  
      if (sum(matches)==0) {
       errorm <- paste("Covar \"",namesaddreg[mnr],"\" in \"addreg.mce\"-object nr ",mcenr,sep="")
       errorm <- paste(errorm,"not in \"covardata\"")
       stop(errorm)
     } else if (sum(matches)==1) {
       mcovar[mnr] <- (1:nrcovardata)[matches]
     } else {
        stop("Dublicated covarnames in \"covardata\"")
     } 
    } 
    
    # We calculate increments
    addregincre <- addregincre[includeincre,]
    machedcovar <- c(as.matrix(covardata[1,mcovar]))
    increhaz[newincrenr,2] <- apply(t(t(addregincre)*machedcovar),1,sum)
    if (estcovar) {
      if (is.vector(addregincrev)) {
        addregincrev <- as.matrix(addregincrev[includeincre])    
        for (obsnr in 1:nrincre) 
          increhaz[newincrenr[obsnr],3] <- t(machedcovar) %*% addregincrev[obsnr] %*% machedcovar   
      } else {
        addregincrev <- addregincrev[includeincre,,]
        for (obsnr in 1:nrincre) 
          increhaz[newincrenr[obsnr],3] <- t(machedcovar) %*% addregincrev[obsnr,,] %*% machedcovar   
      }
    }
  }
   
  
  # We filter transiton matrix data & sort data:
  includetrans <- increhaz[,1]<=min(objinfo[,4])
  if (!is.null(fromtime)) 
    includetrans <- includetrans & (increhaz[,1]>fromtime)
  nrtrans <- sum(includetrans)
  increhaz <- increhaz[includetrans,] 
  increhaz <- increhaz[order(increhaz[,1]),]
  
  # We calculate transition probabilities:
  imat <- diag(maxnrstate)
  transprop <- array(NA,dim=c(nrtrans+1,maxnrstate,maxnrstate))
  pmatrixs  <- array(NA,dim=c(nrtrans+1,maxnrstate,maxnrstate))
  transprop[1,,] <- pmatrixs[1,,] <- imat    
  nrneg <- 0
  nrgt  <- 0
  for (i in 1:nrtrans) {
    pmat <- imat
    from <- increhaz[i,4]
    to   <- increhaz[i,5]
    increh <- increhaz[i,2]
    if (increh<0) {
      nrneg  <- nrneg+1
    }
    if (increh>1) {
      increh <- 1
      nrgt   <- nrgt+1
    }
    pmat[from,to] <- increh
    pmat[from,from] <- 1-increh
    transprop[i+1,,] <- transprop[i,,]%*%pmat
    # We store pmat for use in covariance calculation:
    pmatrixs[i+1,,] <- pmat
  } 
  
  # We calculate covariace:
  if (estcovar) {
    transpvar      <- array(NA,dim=c(nrtrans+1,maxnrstate,maxnrstate))
    transpvar[1,,] <- matrix(0,maxnrstate,maxnrstate)
    
    transpropinv  <- array(NA,dim=c(nrtrans+1,maxnrstate,maxnrstate))
    for (i in 2:(nrtrans+1)) {
      # We find p(u,t) matrix
      transpropinv[i,,] <- imat
      thism <- imat  
      for (j in (i-1):1) {
        thism <- pmatrixs[j+1,,]  %*% thism      
        transpropinv[j,,] <- thism
      }
      nullmat <- diag(rep(0,maxnrstate^2))
      sumcov  <-  matrix(0,maxnrstate^2,maxnrstate^2)
      for (j in 2:i) {
        part2 <- nullmat
        from  <- increhaz[j-1,4]
        to    <- increhaz[j-1,5]
        pnr   <- from+(to-1)*maxnrstate
        pnrf  <- from+(from-1)*maxnrstate
        cvar  <- increhaz[j-1,3]
        part2[pnr,pnr]   <- cvar
        part2[pnrf,pnrf] <- cvar
        part2[pnrf,pnr]  <- (-1)*cvar
        part2[pnr,pnrf]  <- (-1)*cvar        
        part1 <- kronecker(t(transpropinv[j,,]),transprop[j,,])
        part3 <- kronecker(transpropinv[j,,],t(transprop[j,,]))
        sumcov <- sumcov + (part1 %*% part2) %*% part3
      }
      transpvar[i,,] <- matrix(diag(sumcov),maxnrstate)
    } 
  } else {
    transpvar <- NULL
  }

  # We print warnings if needed:
  if (nrneg>0) {
    cat("Remark:",nrneg,"transition estimates gave negative vaules.\n")
  }
  if (nrgt>0) {
    warn <- paste(nrgt,"transitions are truncated because of vaules greater than one")
    warning(warn)
  }

  # We make object and return the result:
  times <- c(increhaz[1,1],increhaz[,1])
  if (is.null(fromtime))
    starttime <- times[1] 
  else 
    starttime <- fromtime
  res <- list(validtrans=vtrans,times=times,transprop=transprop,transpvar=transpvar,
              starttime=starttime,nrnegest=nrneg,nrgt1est=nrgt)
  attr(res,"Call") <- sys.call()
  class(res) <- "addregmc"
  return(res)  
}  

addregmce <- function(addreg.object,from,to) {
  # S-PLUS & R functions for analyses of markov surrviel analyses data with an additive hazard model.
  # Written by:	Harald Fekjaer <harald.fekjar@basalmed.uio.no>, Odd O. Aalen and Ørnulf Borgan
  # See also:	http://www.med.uio.no/imb/stat/addreg/
  # License:	GNU version 2
  # Version:	1.0 Beta 1, © 2001
  #-----------Beginning of routine: -----------------
  
  # We check input data:
  if (!inherits(addreg.object, 'addreg')) 
    stop("Variable \"addreg.object\" must be an addreg object.")
  
  if (from<1 || as.integer(from)!=from)
    stop("Variable \"from\" must be an positive integer.") 

  if (to<1 || as.integer(to)!=to)
    stop("Variable \"to\" must be an positive integer.") 

  if (to==from)
    stop("\"from\" and \"to\" can not be same state.")
   
  if (is.null(addreg.object$increcovar)) {
    errorm <- paste("\"addreg.object\" must include covariance information.")
    errorm <- paste(errorm,"\nPlease run addreg again with \"estcovar=T\"")
    stop(errorm)
  }
    
  # We make object and return the result:
  res <- list(addreg.object=addreg.object,from=from,to=to)
  attr(res,"Call") <- sys.call()
  class(res) <- "addregmce"
  return(res)  
}  

plot.addregmc <- function(addregmc.object,fromto=NULL,stoptime=NULL,confint=F,transtitle=NULL,
                                          xlab="Time",ylab="Probability",standylim=F) {
  # S-PLUS & R functions for analyses of markov surrviel analyses data with an additive hazard model.
  # Written by:	Harald Fekjaer <harald.fekjar@basalmed.uio.no>, Odd O. Aalen and Ørnulf Borgan
  # See also:	http://www.med.uio.no/imb/stat/addreg/
  # License:	GNU version 2
  # Version:	1.0 Beta 1, © 2001
  #-----------Beginning of routine: -----------------
  
  # We check input data:
  if (!inherits(addregmc.object, 'addregmc')) 
    stop("Variable \"addregmc.object\" must be an addregmc object.")
  
  # We read general data:
  validtrans <-addregmc.object$validtrans
  nrtrans <- dim(validtrans)[1]
  if (is.null(stoptime)) 
    lastobs <- length(addregmc.object$times)
  else 
    lastobs <- sum(addregmc.object$times<=stoptime)
  x <- addregmc.object$times[1:lastobs]
    
  # We plot one or several transition propabilits
  if(is.null(fromto)) {
    # We plot several transition propabilits
    # More than one plot are done with recursions 
    if (!is.null(transtitle)) {
      if(length(transtitle)!=nrtrans^2)
        stop("\"transtitle\" must (number of states)^2 number of elements") 
    }
    dimplot <- c(nrtrans,nrtrans)
    save.par <- par("mfrow")
    on.exit(par(mfrow = save.par))
    par(mfrow = dimplot)
    for (tfrom in 1:nrtrans) {
      for (tto in 1:nrtrans) {
        if(is.null(transtitle)) 
          ttitle <- paste(tfrom,"->",tto)
        else
          ttitle <- transtitle[(tfrom-1)*nrtrans+tto]
        if (validtrans[tfrom,tto]) {
           plot.addregmc(addregmc.object,fromto=c(tfrom,tto),stoptime=stoptime,confint=confint,transtitle=ttitle,
                                        xlab=xlab,ylab=ylab,standylim=T)
        } else {
          plot(c(x[1],x[lastobs]),c(0,1),type="n",xlab=xlab,ylab=ylab)
          title(ttitle)
        }
      }     
    }
  } else {
    # We plot one transition propability
    # We chech imput data & start reading data:
    if (!is.vector(fromto))      
      stop("Varible \"fromto\" must be an vector with length 2.")
    if (length(fromto)!=2)
      stop("Varible \"fromto\" must be an vector with length 2")
    tfrom <- fromto[1]
    tto   <- fromto[2]
    if (tfrom>nrtrans || tto>nrtrans)
      stop("\"fromto\" out of transition matrix. (E.g. to large value)")
    if (!validtrans[tfrom,tto])
      stop("\"fromto\" spesifics a not valid transition")
      
    # We read data:
    y <- addregmc.object$transprop[1:lastobs,tfrom,tto]
    if (confint) {
      if (is.null(addregmc.object$transpvar)) {
        errorm <- paste("\"addregmc.object\" must include variance information")
        errorm <- paste(errorm,"\nPlease run \"addregmc\" again with \"estcovar=T\"")
        stop(errorm)
      }
      yse   <- sqrt(addregmc.object$transpvar[1:lastobs,tfrom,tto])
      cintc <- qnorm(.975)
      conf1 <- pmax(0,y-cintc*yse)
      conf2 <- pmin(1,y+cintc*yse)
    } 
        
    # We plot data:
    if (standylim)
      maxintval <- 1
    else if (confint) {
      maxintval <- max(y,conf2)
    } else {
      maxintval <- max(y)
    }
    plot(x,y,type="s",xlab=xlab,ylab=ylab,ylim=c(0,maxintval))
    if (confint) {
      lines(x,conf1,lty = 3,type="s")
      lines(x,conf2,lty = 3,type="s")
    }
    if (is.null(transtitle)) 
       transtitle <- paste("Transitions propabilities from state",tfrom,"to",tto)
     title(transtitle)
  }
}  

plot.addregmce <- function(addregmce.object,...) {
  # S-PLUS & R functions for analyses of markov surrviel analyses data with an additive hazard model.
  # Written by:	Harald Fekjaer <harald.fekjar@basalmed.uio.no>, Odd O. Aalen and Ørnulf Borgan
  # See also:	http://www.med.uio.no/imb/stat/addreg/
  # License:	GNU version 2
  # Version:	1.0 Beta 1, © 2001
  #-----------Beginning of routine: -----------------
  
  # We check input data:
  if (!inherits(addregmce.object, 'addregmce')) 
    stop("Variable \"addregmce.object\" must be an addreg object.")
  
  # We plot the "addreg" part of the object:  
  addreg.object <- addregmce.object$addreg.object
  plot.addreg(addreg.object,...)
}  

"print.addregmc" <- function(addregmc.object) {
  # S-PLUS & R functions for analyses of markov surrviel analyses data with an additive hazard model.
  # Written by:	Harald Fekjaer <harald.fekjar@basalmed.uio.no>, Odd O. Aalen and Ørnulf Borgan
  # See also:	http://www.med.uio.no/imb/stat/addreg/
  # License:	GNU version 2
  # Version:	1.0 Beta 1, © 2001
  #-----------Beginning of routine: -----------------
  
  # We check input data:
  if (!inherits(addregmc.object, 'addregmc')) 
    stop("Variable \"addregmc.object\" must be an addregmc object.")
  
  # We print information about object:  
  cat("Additive hazard regression markov chain object.\n")
  cat("(With markov chain transition propabilites calculated useing Aalen's model)\n\n")
  
  # We print transition matrix:
  vtrans <- addregmc.object$validtrans
  dvtrans <- dim(vtrans)[1]
  vtransprint <- matrix(NA,dvtrans,dvtrans)
  vtransprint[vtrans] <- "X"
  vtransprint[!vtrans] <- " "
  vtransprint <- data.frame(vtransprint)
  names(vtransprint) <- row.names(vtransprint) <- 1:dvtrans
  cat("Valied transitions:\n")
  print(vtransprint)
  cat("\n  Call:\n  ")
  dput(attr(addregmc.object, "Call"))
  invisible(addregmc.object)  
}

print.addregmce <- function(addregmce.object) {
  # S-PLUS & R functions for analyses of markov surrviel analyses data with an additive hazard model.
  # Written by:	Harald Fekjaer <harald.fekjar@basalmed.uio.no>, Odd O. Aalen and Ørnulf Borgan
  # See also:	http://www.med.uio.no/imb/stat/addreg/
  # License:	GNU version 2
  # Version:	1.0 Beta 1, © 2001
  #-----------Beginning of routine: -----------------
  
  # We check input data:
  if (!inherits(addregmce.object, 'addregmce')) 
    stop("Variable \"addregmce.object\" must be an addreg object.")
  
  # We print information about object:  
  addreg.object <- addregmce.object$addreg.object
  cat("Additive hazard regression markov element.\n")
  cat("(With data on one transition useing Aalen's model)\n\n")
  cat("  Transition from:",addregmce.object$from,"\n")
  cat("               to:",addregmce.object$to,"\n\n")
  cat("  Variables:     ")
  cat(dimnames(addreg.object$increments)[[2]])
  cat("\n  Fitted at ")
  cat(addreg.object$ntime[1])
  nint <- length(addreg.object$times)
  cat(" time points until time", addreg.object$times[nint],"\n\n")
  cat("  Call: ")
  dput(attr(addregmce.object, "Call"))
  cat("\n")
  invisible(addregmce.object)  
}  


summary.addregmc <- function(addregmc.object,esttime=NULL,confint=F,digits=3) {
  # S-PLUS & R functions for analyses of markov surrviel analyses data with an additive hazard model.
  # Written by:	Harald Fekjaer <harald.fekjar@basalmed.uio.no>, Odd O. Aalen and Ørnulf Borgan
  # See also:	http://www.med.uio.no/imb/stat/addreg/
  # License:	GNU version 2
  # Version:	1.0 Beta 1, © 2001
  #-----------Beginning of routine: -----------------
  
  # We check input data:
  if (!inherits(addregmc.object, 'addregmc')) 
    stop("Variable \"addregmc.object\" must be an addregmc object.")
  
  cat("\nResultes of markov chain estimation with Additive hazard regression.\n")
  cat("(Estimation of transition propabilites with Aalen's model)\n\n")
  
  # We read data and print transition matrix:
  if (is.null(esttime)) {
    lastobs <- length(addregmc.object$times)
    stopt   <- addregmc.object$times[lastobs]
  } else {
    lastobs <- sum(addregmc.object$times<=esttime)
    stopt   <- esttime
  }
  startt  <- addregmc.object$starttime
  transp  <- addregmc.object$transprop[lastobs,,]
  nrtrans <- dim(transp)[1]
  dimnames(transp) <- list(1:nrtrans,1:nrtrans)
  cat("Transitions propabilities from time",startt,"to",stopt,":\n")
  transp <- round(transp,digits)
  print(transp)
  
  # We print confidens intervals:
  if (confint){
    if (is.null(addregmc.object$transpvar)) {
      errorm <- paste("\"addregmc.object\" must include variance information")
      errorm <- paste(errorm,"\nPlease run \"addregmc\" again with \"estcovar=T\"")
      stop(errorm)
    }
    validtrans <-addregmc.object$validtrans
    nrvtrans   <- sum(validtrans)
    resmat     <- matrix(NA,nrow=nrvtrans,ncol=7)
    resmat[,2] <- rep("->",times=nrtrans)
    vtransnr <- 1  
    transpv  <- addregmc.object$transpvar[lastobs,,]

    # We print transition matrix: 
    for (yfrom in 1:nrtrans) {
      for (xto in 1:nrtrans) {
        if (validtrans[yfrom,xto]) {
          resmat[vtransnr,1] <- yfrom
          resmat[vtransnr,3] <- xto
          obs   <- transp[yfrom,xto]
          obsse <- sqrt(transpv[yfrom,xto])
          cintc <- qnorm(.975)
          resmat[vtransnr,4] <- round(obs,digits)
          resmat[vtransnr,5] <- round(obsse,digits)
          resmat[vtransnr,6] <- max(0,round(obs-cintc*obsse,digits))
          resmat[vtransnr,7] <- min(1,round(obs+cintc*obsse,digits))
          vtransnr <- vtransnr+ 1
        }
      }
    }
    cat("\n")
    resmat <- data.frame(resmat)
    names(resmat) <- c("From"," ","To","Probability","Std. Error","95%","C.I.")
    print(resmat)
  }
  
  # We print call:
  cat("\n  Call:\n  ")
  dput(attr(addregmc.object, "Call"))
}  


summary.addregmce <- function(addregmce.object,...) {
  # S-PLUS & R functions for analyses of markov surrviel analyses data with an additive hazard model.
  # Written by:	Harald Fekjaer <harald.fekjar@basalmed.uio.no>, Odd O. Aalen and Ørnulf Borgan
  # See also:	http://www.med.uio.no/imb/stat/addreg/
  # License:	GNU version 2
  # Version:	1.0 Beta 1, © 2001
  #-----------Beginning of routine: -----------------
  
  # We check input data:
  if (!inherits(addregmce.object, 'addregmce')) 
    stop("Variable \"addregmce.object\" must be an addreg object.")
  
  # We print information about object:  
  addreg.object <- addregmce.object$addreg.object
  cat("Additive hazard regression markov element.\n")
  cat("(With data on one transition useing Aalen's model)\n\n")
  cat("  Transition from:",addregmce.object$from,"\n")
  cat("               to:",addregmce.object$to,"\n\n")
  summary.addreg(addreg.object,...)
  cat("  Call: ")
  dput(attr(addregmce.object, "Call"))
  cat("\n")
  invisible(addregmce.object)  
}