R-scripts/cdd.R

   1 ## PsN R script for plotting results of cdd
   2 ## Justin Wilkins
   3 ## November 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 min.failed    <- FALSE     # do we want to keep minimization failed runs?
  10 cov.failed    <- FALSE      # do we want to keep covariance failed runs?
  11 cov.warnings  <- TRUE      # do we want to keep covariance warning runs?
  12 boundary      <- TRUE      # do we want to keep boundary runs?
  13
  14 ## autogenerated code ends
  15
  16 ## read files
  17 cdd.data <- read.csv("raw_results1.csv")
  18 cdd.inds <- read.csv("skipped_individuals1.csv", header=F)
  19
  20 #cdd.data    <- subset(cdd.data, !is.na(ofv))
  21 #cdd.skipped <- subset(cdd.data, is.na(ofv))
  22
  23 names(cdd.inds)[1] <- "ID"
  24
  25 ## skip first record (no cases deleted)
  26 cdd.data <- cdd.data[-1,]
  27
  28 cdd.data <- cbind(cdd.data,cdd.inds[1])
  29
  30 ## create data frame for current parameter
  31 p1 <- cdd.data
  32
  33 names(p1)[1] <- "minimization.successful"
  34 names(p1)[2] <- "covariance.step.successful"
  35 names(p1)[3] <- "covariance.step.warnings"
  36 names(p1)[4] <- "estimate.near.boundary"
  37
  38 if (min.failed) {
  39   mf  <- subset(p1, minimization.successful == 0)
  40   p1  <- subset(p1, minimization.successful == 1)
  41 }
  42 if (cov.failed) {
  43   cf <- subset(p1, covariance.step.successful == 0)
  44   p1 <- subset(p1, covariance.step.successful == 1)
  45 }
  46 if (!cov.warnings) {
  47   cw <- subset(p1, covariance.step.warnings == 1)
  48   p1 <- subset(p1, covariance.step.warnings == 0)
  49 }
  50 if (!boundary) {
  51   nb <- subset(p1, estimate.near.boundary == 1)
  52   p1 <- subset(p1, estimate.near.boundary == 0)
  53 }
  54
  55 if (exists("mf")) {
  56   cdd.warn <- mf
  57 }
  58 if (exists("cf")) {
  59   if (exists("cdd.warn")) {
  60     cdd.warn <- data.frame(cdd.warn, cf)
  61   } else {
  62     cdd.warn <- cf
  63   }
  64 }
  65 if (exists("cw")) {
  66   if (exists("cdd.warn")) {
  67     cdd.warn <- data.frame(cdd.warn, cw)
  68   } else {
  69     cdd.warn <- cw
  70   }
  71 }
  72 if (exists("nb")) {
  73   if (exists("cdd.warn")) {
  74     cdd.warn <- data.frame(cdd.warn, nb)
  75   } else {
  76     cdd.warn <- nb
  77   }
  78 }
  79
  80 cdd.txt  <- subset(p1, outside.n.sd == 1)
  81 cdd.pt   <- subset(p1, outside.n.sd == 0)
  82
  83 ## Construct plot
  84
  85 postscript(file="cdd.ps", paper="a4",
  86   title="Case-deletion diagnostics",
  87   horizontal=T)
  88
  89 plot (cdd.pt$cov.ratios, cdd.pt$cook.scores,
  90   type="p",
  91   xlab="Cook score",
  92   ylab="Covariance ratio",
  93   main="Case-deletion diagnostics",
  94   xlim=c(0,max(cdd.data$cook.scores, na.rm=T)),
  95   ylim=c(0,max(cdd.data$cov.ratio, na.rm=T))
  96 )
  97
  98 text(cdd.txt$cook.scores, cdd.txt$cov.ratios, labels=as.character(cdd.txt$ID),cex=.8, col=2)
  99 if (exists("cdd.warn")) {
 100   text(cdd.warn$cook.scores, cdd.warn$cov.ratios, labels=as.character(cdd.warn$ID),cex=.8, col=4)
 101 }
 102
 103