R-scripts/llp.R

   1 ## PsN R script for plotting results of llp
   2 ## Justin Wilkins
   3 ## October 2005
   4
   5 ## set basic options
   6 ## autogenerated code begins
   7 ## we can set other things, like colours, line types, layout, etc here if we wish
   8
   9 ofv       <- 3.843     # set this to the desired llp ofv tolerance
  10 ci.lines  <- TRUE      # do we want grey broken lines for outer CI limits?
  11 dofv.line <- TRUE      # do we want grey broken line for final model parameter estimate?
  12 hline     <- TRUE      # do we want grey broken line for desired llp ofv tolerance?
  13
  14 ## autogenerated code ends
  15
  16 ## read file
  17 llp.data <- read.csv("llplog1.csv")
  18
  19 ## get number of parameters
  20 n <- (length(colnames(llp.data)) - 1) / 2
  21
  22 pos <- 0
  23 cols <- c()
  24
  25
  26 ## plots
  27 for (i in 1:(n*2)) {
  28   if (pos != 1) {
  29
  30     pos <- 1
  31     ## create data frame for each parameter
  32     p1 <- c(llp.data[i],llp.data[i+1])
  33     p1 <- as.data.frame(p1)
  34     p1 <- p1[!(apply(is.na(p1),1,any)),]
  35
  36     parmlabel <- names(p1)[1]
  37
  38     names(p1) <- c("Parameter", "dOFV")
  39
  40     ## 3rd order polynomial
  41     fit  <- lm(p1$dOFV~p1$Parameter+I(p1$Parameter^2)+I(p1$Parameter^3))
  42     ## 4th order polynomial
  43     fit2 <- lm(p1$dOFV~p1$Parameter+I(p1$Parameter^2)+I(p1$Parameter^3)+I(p1$Parameter^4))
  44
  45     if (max(p1$Parameter) > 0) {
  46       pmn <- min(p1$Parameter) - (max(p1$Parameter)/10)
  47       pmx <- max(p1$Parameter) + (max(p1$Parameter)/10)
  48     } else {
  49       pmn <- min(p1$Parameter) + (max(p1$Parameter)/10)
  50       pmx <- max(p1$Parameter) - (max(p1$Parameter)/10)
  51     }
  52
  53     pno <- 1000
  54     pdif <- pmx - pmn
  55     pit  <- pdif / pno
  56     pcur <- pmn
  57     bint1 <- FALSE
  58     bint2 <- FALSE
  59
  60     pint1 <- NULL
  61     pint2 <- NULL
  62
  63     ip <- 0 ## inflections - down is bad
  64     misfit <- 0
  65
  66     pfn <- data.frame()
  67     mat.txt <- c()
  68
  69     ## try third order first
  70     for (p in 1:pno) {
  71       pfit <- as.vector( fit$coefficients[1] +
  72           (fit$coefficients[2] * pcur) +
  73           (fit$coefficients[3] * (pcur^2)) +
  74           (fit$coefficients[4] * (pcur^3))
  75         )
  76       addrow  <- c(pcur, pfit)
  77       mat.txt <- c(mat.txt, addrow)
  78
  79       if (p == 1) {
  80         pint1 <- pfit
  81         pint2 <- pfit
  82         pprev <- pfit
  83         tprev <- pcur
  84       } else {
  85         ## going down
  86         if ((pfit <= ofv) && (pprev >= ofv)) {
  87           if (abs(pfit-ofv) <= abs(pprev-ofv)) {
  88             if (!bint1) {
  89               pint1 <- pcur
  90             }
  91           } else {
  92             pint1 <- tprev
  93             ## we have a number for int1
  94             bint1 <- TRUE
  95           }
  96         }
  97         ## going up
  98         if ((pfit >= ofv) && (pprev <= ofv)) {
  99           if (abs(pfit-ofv) <= abs(pprev-ofv)) {
 100             if (!bint2) {
 101               pint2 <- pcur
 102             }
 103           } else {
 104             pint2 <- tprev
 105             ## we have a number for int2
 106             bint2 <- TRUE
 107           }
 108         }
 109         ## have we reached the nadir of the parabola?
 110         if (pfit >= pprev) {
 111           ip <- 1
 112         }
 113         if ((ip == 1) && (pfit < pprev)) {
 114           ## misfit
 115           misfit <- 1
 116         }
 117         pprev <- pfit
 118         tprev <- pcur
 119       }
 120       pcur <- pcur + pit
 121     }
 122
 123     if (misfit) {
 124       ## try fourth order fit
 125       pno <- 1000
 126       pdif <- pmx - pmn
 127       pit  <- pdif / pno
 128       pcur <- pmn
 129       bint1 <- FALSE
 130       bint2 <- FALSE
 131
 132       pint1 <- NULL
 133       pint2 <- NULL
 134
 135       pfn <- data.frame()
 136       mat.txt <- c()
 137
 138       for (p in 1:pno) {
 139         pfit <- as.vector( fit2$coefficients[1] +
 140             (fit2$coefficients[2] * pcur) +
 141             (fit2$coefficients[3] * (pcur^2)) +
 142             (fit2$coefficients[4] * (pcur^3)) +
 143             (fit2$coefficients[5] * (pcur^4))
 144           )
 145         addrow  <- c(pcur, pfit)
 146         mat.txt <- c(mat.txt, addrow)
 147
 148         if (p == 1) {
 149           pint1 <- pfit
 150           pint2 <- pfit
 151           pprev <- pfit
 152           tprev <- pcur
 153         } else {
 154           ## going down
 155           if ((pfit <= ofv) && (pprev >= ofv)) {
 156             if (abs(pfit-ofv) <= abs(pprev-ofv)) {
 157               if (!bint1) {
 158                 pint1 <- pcur
 159               }
 160             } else {
 161               pint1 <- tprev
 162               ## we have a number for int1
 163               bint1 <- TRUE
 164             }
 165           }
 166           ## going up
 167           if ((pfit >= ofv) && (pprev <= ofv)) {
 168             if (abs(pfit-ofv) <= abs(pprev-ofv)) {
 169               if (!bint2) {
 170                 pint2 <- pcur
 171               }
 172             } else {
 173               pint2 <- tprev
 174               ## we have a number for int2
 175               bint2 <- TRUE
 176             }
 177           }
 178           ## have we reached the nadir of the parabola?
 179           if (pfit >= pprev) {
 180             ip <- 1
 181           }
 182           if ((ip == 1) && (pfit < pprev)) {
 183             ## misfit
 184             misfit <- 1
 185           }
 186           pprev <- pfit
 187           tprev <- pcur
 188         }
 189         pcur <- pcur + pit
 190       }
 191     }
 192
 193     mdat <- matrix(mat.txt, nrow = pno, ncol=2, byrow=TRUE,
 194       dimnames = list(c(), c("Parameter", "dOFV")))
 195
 196     llp.parabola <- as.data.frame(mdat)
 197
 198     postscript(file=paste("llp.", parmlabel, ".ps", sep=""), paper="a4",
 199       title=paste("Log-likelihood profile - ", parmlabel, sep=""),
 200       horizontal=T)
 201
 202     plot (llp.parabola$Parameter, llp.parabola$dOFV,
 203       type="l",
 204       # ylim=c(0,max(data.main$DV)),
 205       ylab="dOFV",
 206       xlab=parmlabel,
 207       main=paste("Log-likelihood profile - ", parmlabel, sep="")
 208       #sub="Binning tolerance: 0.25, occasion: 0"
 209       )
 210
 211     points(p1$Parameter, p1$dOFV, col="red", pch=21)
 212     abline(v=pint1, col="red", lty=2)
 213     abline(v=pint2, col="red", lty=2)
 214     text(pint1, ofv, labels=signif(pint1, digits = 3), cex = .8, adj = c(0,0), pos='4')
 215     text(pint2, ofv, labels=signif(pint2, digits = 3), cex = .8, adj = c(0,0), pos='4')
 216     text(p1[1,1], 0, labels=signif(p1[1,1], digits = 3), cex = .8, adj = c(0,0), pos='3')
 217
 218     if (ci.lines) {
 219       abline(v=min(p1$Parameter), col="gray", lty=2)
 220       abline(v=max(p1$Parameter), col="gray", lty=2)
 221     }
 222
 223     if (dofv.line) {
 224       abline(v=p1[1,1], col="gray", lty=2)
 225     }
 226
 227     if (hline) {
 228       abline(h=ofv, col="gray", lty=2)
 229     }
 230
 231   } else {
 232     pos <- 0
 233   }
 234 }
 235
 236