R/reportHTML.spikeTrain.R

In STAR: Spike Train Analysis with R

reportHTML.spikeTrain <- function(object,
                                  filename,
                                  extension="html",
                                  directory=getwd(),
                                  Title,
                                  forceTT=TRUE,
                                  digits=3,
                                  timeUnit="s",
                                  otherST,
                                  laglim=c(-0.1,0.1),
                                  cch=c("both","scch","cch"),
                                  doGamCheck=TRUE,
                                  k=100,
                                  bs="tp",
                                  nbEvtPerBin=50,
                                  ...
                                  )
### reportHTML method for spikeTrain objects.
### A spike train plot of object is generated first
### by a call to plot.spikeTrain. A short numerical
### summary is generated next (partial output of a call
### to summary.spikeTrain). A renewal test is performed
### by a call to renewalTestPlot using default arguments
### of the latter. The six duration distributions with
### two paramters are fitted with a call to compModels with
### the plot argument set to FALSE and a call to xMLE where
### "x" stands for the best model (in the AIC sense) among the
### the 6. A summary (best value and se of the two model 
### parameters) of the best fit is printed. A time transformation using:
### cumsum(px(diff(y))
### where "x" as the same meaning as above and y stands for the
### time transformed (TT) version of the original train using the
### best among the 6 models. If this TT version passes the "tests
### of Ogata" at the 99% level (these tests are returned by a call
### to summary.transformedTrain), then a plot is generated by a
### new call to compModels displaying the QQ plot on a log scale
### for each of the 6 fitted models as well as the Ogata's test
### plots. If arg forceTT is set to TRUE, these plots are generated
### even if the tests are not passed.
### Args "digits" and "timeUnit" are the corresponding args of
### summary.spikeTrain
{

  objectN <- deparse(substitute(object))
  
  ## check is object is a spikeTrain object
  if (!is.spikeTrain(object)) object <- as.spikeTrain(object)

  if (missing(filename))
    filename <- paste(objectN,"analysis")
  
  if (missing(Title))
    Title <- filename

  
  HTMLInitFile(outdir=directory,
               filename=filename,
               extension=extension,
               Title=Title)

  fullName <- paste(directory,"/",
                    filename,".",
                    extension,sep="")

  saveName <- paste(directory,"/",
                    filename,".rda",
                    sep="")
  
  HTML.title(filename,
             file=fullName,
             HR=2)
  
  ## add a spike train plot 
  HTML.title(paste("Spike train plot of ",
                   filename,sep=""),
             file=fullName,HR=3)

  
  stFigName <- paste(filename,"_st.png",sep="")
  figFname <- paste(directory,"/",stFigName,sep="")
  png(figFname,width=500,height=500)
  plot(object)
  dev.off()
  HTMLInsertGraph(stFigName,
                  file=fullName,
                  WidthHTML=500,
                  HeightHTML=500)

  HTMLbr(2,file=fullName)
  HTML.title(paste("Short summary of ",
                   filename,sep=""),
             file=fullName,HR=3)
  stRange <- range(object)
  stNb <- length(object)
  isi <- diff(object)
  stStat1 <- c(mean(isi), sd(isi))
  stStat2 <- c(mean(log(isi)), sd(log(isi)))
  cat(paste("<p>A spike train with ", stNb, " events, starting at: ", 
            round(stRange[1], digits = digits), " and ending at: ", 
            round(stRange[2], digits = digits), " (", timeUnit, ").</p>", 
            "<p>The mean ISI is: ", round(stStat1[1], digits = digits), 
            " and its SD is: ", round(stStat1[2], digits = digits), 
            " (", timeUnit, ").</p>", "<p>The shortest interval is: ", 
            round(min(isi), digits = digits), " and the longest is: ", 
            round(max(isi), digits = digits), " (", timeUnit, ").</p>", 
            sep = ""),
      file=fullName,
      append=TRUE)
  
  ## add a renewal test plot if there are more than 50 events
  if (length(object) > 50) {
    HTMLhr(file=fullName)
    HTML.title(paste("Renewal test of ",
                     filename,sep=""),
               file=fullName,HR=3)
    
    rtFigName <- paste(filename,"_rt.png",sep="")
    figFname <- paste(directory,"/",rtFigName,sep="")
    png(figFname,width=800,height=800)
    renewalTestPlot(object)
    dev.off()
    HTMLInsertGraph(rtFigName,
                    file=fullName,
                    WidthHTML=800,
                    HeightHTML=800)
  } ## End of conditional on length(object) > 50

  ## compare models
  cm <- compModels(object,plot=FALSE)
  ## get best model
  bestM <- names(cm)[1]
  HTMLhr(file=fullName)
  cat(paste("<p>The best model with 2 parameters is: ",
            bestM,
            ".</p>",sep=""),
      file=fullName,
      append=TRUE
      )
  cat("<p>Its estimated parameters values and associated se are:</p>",
      file=fullName,
      append=TRUE
      )
  
  bestFit <- switch(bestM,
                    invgauss=invgaussMLE(object),
                    lnorm=lnormMLE(object),
                    gamma=gammaMLE(object),
                    weibull=weibullMLE(object),
                    llogis=llogisMLE(object),
                    rexp=rexpMLE(object)
                    )
  
  toPrint <- rbind(bestFit$estimate,bestFit$se)
  rownames(toPrint) <- c("mean","se")
  HTML(toPrint,file=fullName)
  
  Lambda <- switch(bestM,
                   invgauss=pinvgauss(isi,bestFit$estimate[1],
                     bestFit$estimate[2],log.p=TRUE,lower.tail=FALSE),
                   lnorm=plnorm(isi,bestFit$estimate[1],
                     bestFit$estimate[2],log.p=TRUE,lower.tail=FALSE),
                   gamma=pgamma(isi,bestFit$estimate[1],
                     bestFit$estimate[2],log.p=TRUE,lower.tail=FALSE),
                   weibull=pweibull(isi,bestFit$estimate[1],
                     bestFit$estimate[2],log.p=TRUE,lower.tail=FALSE),
                   llogis=pllogis(isi,bestFit$estimate[1],
                     bestFit$estimate[2],log.p=TRUE,lower.tail=FALSE),
                   rexp=prexp(isi,bestFit$estimate[1],
                     bestFit$estimate[2],log.p=TRUE,lower.tail=FALSE)
                   )
  Lambda <- -cumsum(Lambda)
  class(Lambda) <- c("transformedTrain","spikeTrain")
  
  if (max(Lambda)/50 > 2) {
    LambdaS <- summary(Lambda)

    otTRUE <- LambdaS[[1]][2] &&
    LambdaS[[2]][2] &&
    (0.005 <= LambdaS[[3]][2]) &&
    (0.995 >= LambdaS[[3]][2])
  } else {
    otTRUE <- FALSE
  } ## End of conditional on max(Lambda)/50 > 2
  
  if (otTRUE || forceTT) {
    HTMLbr(1,file=fullName)
    HTML.title(paste("Model comparison for ",
                     filename,sep=""),
               file=fullName,HR=3)
    cmFigName <- paste(filename,"_cm.png",sep="")
    figFname <- paste(directory,"/",cmFigName,sep="")
    png(figFname,width=800,height=800)
    compModels(object)
    dev.off()
    HTMLInsertGraph(cmFigName,
                    file=fullName,
                    WidthHTML=800,
                    HeightHTML=800)

    HTMLbr(3,file=fullName)
    HTML.title(paste("Ogata's tests after time transformation of ",
                     filename,
                     " using a ",bestM," model",
                     sep=""),
             file=fullName,HR=3)
  
    otFigName <- paste(filename,"_TTot.png",sep="")
    figFname <- paste(directory,"/",otFigName,sep="")
    png(figFname,width=800,height=800)
    if (max(Lambda)/50 > 2) {
      plot(Lambda,
           which=c(1,2,4,5),
           ask=FALSE)
    } else {
      plot(Lambda,
           which=c(1,2,4),
           ask=FALSE)
    }
    dev.off()
    HTMLInsertGraph(otFigName,
                    file=fullName,
                    WidthHTML=800,
                    HeightHTML=800)
  } ## End of conditional on otTRUE 

  ## check if otherST is given
  if (!missing(otherST)) {

    keepGoing <- TRUE
    ## check that otherST is a named list of spikeTrain objects
    if (!inherits(otherST,"list") || is.null(names(otherST))) {
      warning("otherST should be a named list of spikeTrain objects.")
      keepGoing <- FALSE
    }
    
    if (keepGoing) {

      if (cch[1] %in% c("both","scch")) {
        scchL <- lapply(otherST,
                        function(st) {
                          lt <- lockedTrain(object,st,laglim=laglim)
                          testF <- lt$nbTestSpikes/lt$obsTime
                          nRef <- lt$nbRefSpikes
                          theBW <- round(nbEvtPerBin/testF/nRef,digits=3)
                          k <- min(k,diff(lt$laglim)%/%theBW)
                          slockedTrain(lt,bw=theBW,bs=bs,k=k)
                        }
                        )
      } ## End of conditional on cch[1] %in% c("both","scch")

      if (cch[1] %in% c("both","cch")) {
        cchL <- lapply(otherST,
                       function(st) {
                         lt <- lockedTrain(object,st,laglim=laglim)
                         testF <- lt$nbTestSpikes/lt$obsTime
                         nRef <- lt$nbRefSpikes
                         theBW <- round(nbEvtPerBin*1.5/testF/nRef,digits=4)
                         hist(lt,bw=theBW,plot=FALSE)
                        }
                        )
      } ## End of conditional on cch[1] %in% c("both","cch")

      sapply(seq(length(otherST)),
             function(trainIdx) {

               ## Write general header
               HTMLbr(1,file=fullName)
               HTMLhr(file=fullName)
               HTML.title(paste("Cross-Intensity with ",
                                names(otherST)[trainIdx],
                                " as a test train",sep=""),
                          file=fullName,HR=3)
               
               if (cch[1] %in% c("both","scch")) {
                 ## A smoothed cross-intensity was evaluated
                 ## write down its summary
                 HTML.title(paste("GAM fit summary ",
                                  "(pre-binning bin width, ",
                                  scchL[[trainIdx]][["bwV"]][1],
                                  ", k=",
                                  length(coef(scchL[[trainIdx]][["gamFit"]]))-1,
                                  ", bs=\"",bs,"\"):",sep=""),
                            file=fullName,HR=4)
                 fitS <- summary(scchL[[trainIdx]])
                 HTML(fitS,file=fullName)

                 if (doGamCheck) {  
                   HTMLbr(1,file=fullName)
                   HTML.title("GAM goodness of fit diagnostics:",,
                              file=fullName,HR=4)
                   gcFigName <- paste(filename,
                                      "_",names(otherST)[trainIdx],
                                      "_gc.png",sep="")
                   figFname <- paste(directory,"/",gcFigName,sep="")
                   png(figFname,width=800,height=800)
                   gam.check(scchL[[trainIdx]][["gamFit"]])
                   dev.off()
                   HTMLInsertGraph(gcFigName,
                                   file=fullName,
                                   WidthHTML=800,
                                   HeightHTML=800)
                 } ## End of conditional on doGamCheck
                 
                 if (cch[1] %in% c("both")) {
                   ## plot both cross-intensity estiamtes
                   HTML.title("Smoothed and \"classical\" cross-intensity plots:",
                              file=fullName,HR=4)
                   ciFigName <- paste(filename,
                                      "_",names(otherST)[trainIdx],
                                      "_ci.png",sep="")
                   figFname <- paste(directory,"/",ciFigName,sep="")
                   png(figFname,width=1000,height=700)
                   layout(matrix(1:2,1,2))
                   plot(scchL[[trainIdx]],main="")
                   plot(cchL[[trainIdx]],main="")
                   dev.off()
                   HTMLInsertGraph(ciFigName,
                                   file=fullName,
                                   WidthHTML=1000,
                                   HeightHTML=700)
                 } else {
                   ## plot only the smoothed intensity estimate
                   HTML.title("Smoothed cross-intensity plot:",
                              file=fullName,HR=4)
                   ciFigName <- paste(filename,
                                      "_",names(otherST)[trainIdx],
                                      "_ci.png",sep="")
                   figFname <- paste(directory,"/",ciFigName,sep="")
                   png(figFname,width=500,height=500)
                   plot(scchL[[trainIdx]],main="")
                   dev.off()
                   HTMLInsertGraph(ciFigName,
                                   file=fullName,
                                   WidthHTML=500,
                                   HeightHTML=500)
                 } ## End of conditional on cch[1] %in% c("both")
               } else {
                 ## plot the classical cross-intensity estimate
                 HTML.title("\"Classical\" cross-intensity plots:",
                            file=fullName,HR=4)
                 ciFigName <- paste(filename,
                                    "_",names(otherST)[trainIdx],
                                    "_ci.png",sep="")
                 figFname <- paste(directory,"/",ciFigName,sep="")
                 png(figFname,width=500,height=500)
                 plot(cchL[[trainIdx]],main="")
                 dev.off()
                 HTMLInsertGraph(ciFigName,
                                 file=fullName,
                                 WidthHTML=500,
                                 HeightHTML=500)
               } ## End of conditional on cch[1] %in% c("both","scch")
             }
             )
      
    } ## End of conditional on keepGoing
    
  } ## End of conditional on !missing(otherST)
  
  HTMLEndFile()

  fctCall <- match.call()

  if (!exists("scchL") && !exists("cchL")) { 
    save(cm,bestFit,Lambda,fctCall,file=saveName)
  } else {
    if (exists("scchL") && exists("cchL")) {
      save(cm,bestFit,Lambda,fctCall,scchL,cchL,file=saveName)
    } else {
      if (exists("scchL")) {
        save(cm,bestFit,Lambda,fctCall,scchL,file=saveName)
      } else {
        save(cm,bestFit,Lambda,fctCall,cchL,file=saveName)
      } ## End of conditional on exists("scchL")
    } ## End of conditional on exists("scchL") && exists("cchL")
  } ## End of conditional on !exists("scchL") && !exists("cchL")
  
}