R-scripts/cdd.R

   1 ## PsN R script for plotting results of cdd
   2 ## Justin Wilkins
   3 ## November 2005
   4 ## Lars Lindbom
   5 ## August 2006
   6
   7
   8 ## set basic options
   9 ## autogenerated code begins
  10 ## we can set other things, like colours, line types, layout, etc here if we wish
  11
  12 min.failed    <- FALSE     # do we want to keep minimization failed runs?
  13 cov.failed    <- FALSE      # do we want to keep covariance failed runs?
  14 cov.warnings  <- TRUE      # do we want to keep covariance warning runs?
  15 boundary      <- TRUE      # do we want to keep boundary runs?
  16
  17 # We have three options for the markers:
  18 # 1. as circles
  19 # 2. as text (taken from the column identified by the -case_column
  20 #    option to the PsN cdd command)
  21 # 3. as circles and text. The subset of runs which fall outside of 2
  22 #    standard deviations around the center of the cook-scores and
  23 #    covariance ratios as calculated using principal components
  24 #    analysis are marked as in 2.
  25
  26 markeropt <- 2
  27
  28 ## autogenerated code ends
  29
  30 ## read files
  31 cdd.data <- read.csv("raw_results1.csv")
  32 cdd.inds <- read.csv("skipped_individuals1.csv", header=F)
  33
  34 #cdd.data    <- subset(cdd.data, !is.na(ofv))
  35 #cdd.skipped <- subset(cdd.data, is.na(ofv))
  36
  37 names(cdd.inds)[1] <- "ID"
  38
  39 ## skip first record (no cases deleted)
  40 cdd.data <- cdd.data[-1,]
  41
  42 cdd.data <- cbind(cdd.data,cdd.inds[1])
  43
  44 ## create data frame for current parameter
  45 p1 <- cdd.data
  46
  47 names(p1)[1] <- "minimization.successful"
  48 names(p1)[2] <- "covariance.step.successful"
  49 names(p1)[3] <- "covariance.step.warnings"
  50 names(p1)[4] <- "estimate.near.boundary"
  51
  52 if (min.failed) {
  53   mf  <- subset(p1, minimization.successful == 0)
  54   p1  <- subset(p1, minimization.successful == 1)
  55 }
  56 if (cov.failed) {
  57   cf <- subset(p1, covariance.step.successful == 0)
  58   p1 <- subset(p1, covariance.step.successful == 1)
  59 }
  60 if (!cov.warnings) {
  61   cw <- subset(p1, covariance.step.warnings == 1)
  62   p1 <- subset(p1, covariance.step.warnings == 0)
  63 }
  64 if (!boundary) {
  65   nb <- subset(p1, estimate.near.boundary == 1)
  66   p1 <- subset(p1, estimate.near.boundary == 0)
  67 }
  68
  69 if (exists("mf")) {
  70   cdd.warn <- mf
  71 }
  72 if (exists("cf")) {
  73   if (exists("cdd.warn")) {
  74     cdd.warn <- data.frame(cdd.warn, cf)
  75   } else {
  76     cdd.warn <- cf
  77   }
  78 }
  79 if (exists("cw")) {
  80   if (exists("cdd.warn")) {
  81     cdd.warn <- data.frame(cdd.warn, cw)
  82   } else {
  83     cdd.warn <- cw
  84   }
  85 }
  86 if (exists("nb")) {
  87   if (exists("cdd.warn")) {
  88     cdd.warn <- data.frame(cdd.warn, nb)
  89   } else {
  90     cdd.warn <- nb
  91   }
  92 }
  93
  94 if( markeropt == 1 ) {
  95   cdd.txt  <- subset(p1, FALSE)
  96   cdd.pt   <- subset(p1, TRUE)
  97 }
  98
  99 if( markeropt == 2 ) {
 100   cdd.txt  <- subset(p1, TRUE)
 101   cdd.pt   <- subset(p1, FALSE)
 102 }
 103
 104 if( markeropt == 3 ) {
 105   cdd.txt  <- subset(p1, outside.n.sd == 1)
 106   cdd.pt   <- subset(p1, outside.n.sd == 0)
 107 }
 108
 109
 110 ## Construct plot
 111
 112 pdf(file="cdd.pdf", paper="special",
 113   title="Case-deletion diagnostics",width=10,height=7)
 114
 115 plot (cdd.pt$cov.ratios, cdd.pt$cook.scores,
 116   type="p",
 117   xlab="Cook score",
 118   ylab="Covariance ratio",
 119   main="Case-deletion diagnostics",
 120   xlim=c(0,max(cdd.data$cook.scores, na.rm=T)),
 121   ylim=c(0,max(cdd.data$cov.ratio, na.rm=T))
 122 )
 123
 124 text(cdd.txt$cook.scores, cdd.txt$cov.ratios, labels=as.character(cdd.txt$ID),cex=.8, col=2)
 125 if (exists("cdd.warn")) {
 126   text(cdd.warn$cook.scores, cdd.warn$cov.ratios, labels=as.character(cdd.warn$ID),cex=.8, col=4)
 127 }
 128
 129