#! /usr/bin/R
sink(file="/dev/null") # to stop getting the Hmisc announcement text
graphDir <- "/var/psyctc/psyctc.org/stats/R/CSC1/scratch/"
graphURLroot <- "/stats/R/CSC1/scratch/"
if(length(system(paste("ls ",graphDir,"*.png",sep=""), intern=TRUE)) > 0 )
   system(paste("rm ",graphDir,"*.png",sep=""))


sink() # switch output back to stdout

tag(HTML)
    tag(HEAD)
        tag(TITLE)
            cat("Criteria for &quot;Clinically Significant Change&quot; (CSC)")
        untag(TITLE)
    untag(HEAD)

    lf(2)

    tag(BODY, bgcolor = "lime")
        lf(2)
        tag(center)
        cat("<H1>Criteria for &quot;Clinically Significant Change&quot; (CSC)</H1>")

        comments("Get the data")
        prob <- 0
        cmean <- as.numeric(scanText(formData$cmean))
        ncmean <- as.numeric(scanText(formData$ncmean))
        csd <- as.numeric(scanText(formData$csd))
        ncsd <- as.numeric(scanText(formData$ncsd))
        smin <- as.numeric(scanText(formData$smin))
        smax <- as.numeric(scanText(formData$smax))
        round <- as.numeric(scanText(formData$round))
        large <- as.numeric(scanText(formData$large))
        comments("Test the input data")
        if (is.na(cmean) || is.na(ncmean) || is.na(csd) || is.na(ncsd)) {
           cat("<H2>You must give all four parameters as numeric values, go back and try again!</H2>")
           prob <- 1
         }
        if ((round < 0) | (round > 6)) {
           cat("<H2>Program only offers rounding to between zero and six decimal places.  Go back and try again!</H2>")
           prob <- 1
         }
        if (prob) {
           cat("<H1>Go back and try again!</H1>")
           cat("CGI script written by Chris Evans using David Firth's excellent GGIwithR package.  ")
           linkto("Contact me if something isn't working right ...", "http://www.psyctc.org/cgi-bin/mailto.pl?webmaster") ; br()
         }
        lf(2)
        comments("end of testing")

        if (!prob) {
           comments("Got into results section")
           if (large) {
              cat("The graphic below may look huge but should be good enough to save to your hard disc and use for publications")
              cat(" or presentations.  Some browsers may have rescaled it to make it fit your screen but saving it should")
              cat(" always give you all the detail of the original.  The text results follow below, after the plot.")
           }
           cat("<TABLE BORDER='1' CELLSPACING='1' CELLPADDING='4' WIDTH='620' BGCOLOR='lime')")
           cat("<TR>")
           cat("<TH BGCOLOR='white' COLSPAN='4'>")
           cat("<B>Your input</B>")
           cat("</TH>")
           cat("</TR>")
           cat("<tr>")
           cat("<td colspan='4'>")
              now <- system('date +%Y-%m-%d,%T', intern=TRUE)
              host <- system("echo $REMOTE_ADDR", intern=TRUE)
              cat("Request from: ",host, " at", now,"<BR>")
           cat("</td>")
       cat("</tr>")
       cat("<tr>")
           cat("<th COLSPAN='2'>Clinical</th><th COLSPAN='2'>Non-clinical</th></th>")
       cat("</tr>")
       cat("<tr>")
           cat("<td>")
              cat("Mean")
           cat("</td>")
           cat("<td>")
              cat(cmean)
           cat("</td>")
           cat("<td>")
              cat("Mean")
           cat("</td>")
           cat("<td>")
              cat(ncmean)
           cat("</td>")
       cat("</tr>")
       cat("<tr>")
           cat("<td>")
              cat("s.d.:")
           cat("</td>")
           cat("<td>")
              cat(csd)
           cat("</td>")
           cat("<td>")
              cat("s.d.:")
           cat("</td>")
           cat("<td>")
              cat(ncsd)
           cat("</td>")
       cat("</tr>")
       if (!is.na(smin)) {
          cat("<tr>")
          cat("<td COLSPAN='3'>")
             cat("Minimum possible score on the measure: ")
          cat("</td>")
          cat("<td COLSPAN='1'>")
             cat(smin)
          cat("</td>")
          cat("</tr>")
       }
       if (!is.na(smax)) {
          cat("<tr>")
          cat("<td COLSPAN='3'>")
             cat("Maximum possible score on the measure: ")
          cat("</td>")
          cat("<td COLSPAN='1'>")
             cat(smax)
          cat("</td>")
          cat("</tr>")
       }
       cat("<tr>")
       cat("<td COLSPAN='4'>")
       cat("<p>")
       cat("Output produced at ", date())
       cat("</p>")
       cat("<tr>")
       cat("<th bgcolor='white' COLSPAN='4'>")
       cat("<b>Results for your parameters</b>")
       cat("</th>")
       cat("</tr>")

       ## plots the Jacobson CSC criteria

       if (cmean < ncmean) {
           comments("reversed = 1")
           reversed <- 1
       } else {
           comments("reversed = 0")
           reversed <- 0
       }
       # range of scores possible on the instrument
       if (is.na(smin)) {
           cmin <- cmean - 3*csd
           ncmin <- ncmean - 3*ncsd
           xmin <- min(cmin,ncmin)
       } else {
           xmin <- smin
       }
       if (is.na(smax)) {
           cmax <- cmean + 5*csd
           ncmax <- ncmean + 5*ncsd
           xmax <- max(cmax,ncmax)
       } else {
           xmax <- smax
       }
       inc <- (xmax - xmin)/100 # used for prettified placement of text

       crit.c <- (cmean*ncsd + ncmean*csd)/(ncsd+csd)
       if (reversed) crit.a <- cmean + 2*csd
          else crit.a <- cmean - 2*csd
       if (reversed) crit.b <- ncmean - 2*ncsd
          else crit.b <- ncmean + 2*ncsd
       crit.a.rnd <- round(crit.a,round)
       crit.b.rnd <- round(crit.b,round)
       crit.c.rnd <- round(crit.c,round)
           
       cat("<tr>")
       cat("<td COLSPAN='3'>")
           cat("Criterion A:")
       cat("</td>")
       cat("<td COLSPAN='1'>")
           cat(crit.a.rnd)
       cat("</td>")
       cat("<tr>")
       cat("<td COLSPAN='3'>")
           cat("Criterion B:")
       cat("</td>")
       cat("<td COLSPAN='1'>")
           cat(crit.b.rnd)
       cat("</td>")
       cat("<tr>")
       cat("<td COLSPAN='3'>")
           cat("Criterion C:")
       cat("</td>")
       cat("<td COLSPAN='1'>")
           cat(crit.c.rnd)
       cat("</td>")
       
       file <- paste("tmp",host,".",now,".png",sep="")
       if (large) {
           webPNG(file=file,height=8,width=10,res=360)
       } else {
           webPNG(file=file,height=8,width=10,res=72)
       }
       ## nothing beyond here needs to be tweaked unless you want different plot formatting
       # get the Gaussian quantiles
       x <- seq(xmin,xmax,inc/100) # possible score range on a fairly tight plotting frequency
       nclin <- dnorm(x,ncmean,ncsd)
       clin <- dnorm(x,cmean,csd)
       # set up the plotting space
       plot(nclin,xlim=c(xmin,xmax),ylim=c(0,1),ylab="",xlab="",type="n")
       # put in vertical reference lines on the limits of the possible scores
       par(lwd=1)
       par(lty=2)
       abline(v=xmin)
       abline(v=xmax)
       # put in vertical reference lines on the two means
       par(adj=.5)
       par(lty=4)
       abline(v=ncmean)
       text(ncmean,.9,"Non-clin. mean")
       abline(v=cmean)
       text(cmean,.9,"Clinical mean")
       par(lty=1)
       # double the thickness of the lines and put in the Gaussian quantiles
       par(lty=1)
       par(lwd=2)
       lines(x,nclin)
       lines(x,clin)
       ## change colour (to red as coded here) for the Jacobson criteria
       par(col=2)
       abline(v=crit.c)
       text(crit.c,.8,"CSC criterion C")
       arrows(cmean,.75,crit.c,.75,col=2)
       arrows(ncmean,.75,crit.c,.75,col=2)
       abline(v=crit.a)
       par(adj=1)
       text(crit.a-inc,.6,"CSC criterion A")
       arrows(cmean,.55,crit.a,.55,col=2)
       abline(v=crit.b)
       par(adj=0)
       text(crit.b+inc,.4,"CSC criterion B")
       arrows(ncmean,.35,crit.b,.35,col=2)
       par(lty=1)
       par(adj=.5)
       par(col=1)
    graphics.off()
    img(src = file) ; br(2)




        cat("</center></b>")
        cat("</td>")
        cat("</tr>")
        cat("<tr>")
        cat("<tr>")
        cat("<th bgcolor='white' COLSPAN='4'>")
        cat("<b>Explanation and advice for CORE system users and others</b>")
        cat("</th>")
        cat("</tr>")
        cat("<tr>")
        cat("<td COLSPAN='4'>")
        cat("The logic of this was described in:<BR>")
        cat("&nbsp;&nbsp;&nbsp;Jacobson, N. S., Follette, W. C. & Revenstorf, D. (1984) Psychotherapy outcome research: methods for reporting variability and evaluating clinical significance. Behavior Therapy, 15, 336-352.<BR>")
        cat("and nicely summarised in:<BR>")
        cat("&nbsp;&nbsp;&nbsp;Jacobson, N. S. & Truax, P. (1991) Clinical significance: a statistical approach to defining meaningful change in psychotherapy research. Journal of Consulting and Clinical Psychology, 59, 12-19.<BR>")
        cat("and even more succinctly in:<BR>")
        cat("&nbsp;&nbsp;&nbsp;Evans, C., Margison, F. & Barkham, M. (1998) The contribution of reliable and clinically significant change methods to evidence-based mental health. Evidence Based Mental Health, 1, 70-72.<BR>")
        cat("<P><UL><LI>Criterion A is that the post-treatment score should like more than two standard deviations from the mean of a clinical population")
        cat("<LI>Criterion B is that the post-treatment score should lie within two standard deviations of the &quot;normal&quot; population mean")
        cat("<LI>Criterion C is that the score should go across a point halfway between the two means.  It's actually halfway in terms of the spreads of the two distributions, it won't be exactly halfway if the spreads of the &quot;normal&quot; and clinical distributions are different.  In that situation it will be an equal distance from the two means <i>in terms of SD units</i>, i.e. the criterion will be nearer the mean with the smaller SD.  This preserves the wonderful property of criterion C, on Gaussian data (aye, there's the rub!) of giving you equal misclassification rates: the proportion of the clinical population classified as normal, false negatives, is equal to the proportion of the normal population misclassified as cases, false positives.")
        cat("</UL>Jacobson and Truax (p.14) suggest that where &quot;normal&quot; distributions are known, and the two populations overlap,")
        cat("criterion C is best but B is better if the two distributions ae non-overlapping and they note that when non-clinical")
        cat("distribution data is not available, A is the only criterion that can be used.")
        cat("<P>All this is based on the model of Gaussian (infinite) population distributions.")
        cat(" The maths are probably robust to finite populations but will break down for")
        cat("markedly non-Gaussian distributions, i.e. the real world.")
        cat("</td>")
        cat("</tr>")

        cat("<tr>")
        cat("<th bgcolor='white' COLSPAN='4'>")
        cat("<b>Technicalities</b>")
        cat("</th>")
        cat("</tr>")
        cat("<tr>")
        cat("<td COLSPAN='4'>")
        cat("Calculation done in R, CGI script written by Chris Evans using David Firth's excellent GGIwithR package which can be obtained from CRAN..  ")
     linkto("Contact me if something isn't working right ...", "http://www.psyctc.org/cgi-bin/mailto.pl?webmaster") ; br()

     lf()

     lf()
     cat("</td>")
     cat("</tr>")
     cat("</table>")
     }
  untag(BODY)
untag(HTML)
lf()
sink(paste("/home/xychris0/R/CGIwithR-log/CSC1.R.",as.character(now),sep=""))
cat("program = CSC1.R ; host = ",host, "; ", now, "; cmean = ",cmean,"; ncmean = ",ncmean,"; csd = ",csd,\
    "; ncsd = ",ncsd,"; round = ",round,"\n")
sink()