R-scripts/llp.R

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