add JF extra/confidenceEllipses.R

extra/Toy.R

@@ -145,6 +145,20 @@ car::confidenceEllipse(Toy.mlm, which=c(2,4), levels = 0.68,
                        fill = TRUE, fill.alpha = 0.2,
                        cex.lab = 1.5)
 
+heplots::coefplot(Toy.mlm)
+# FIXME
+# Error in coefplot.mlm(Toy.mlm) : There are only 0 response variables.
+
+# try JF's new version
+source(here("extra", "confidenceEllipses.R"))
+
+data(Duncan, package="carData")
+m <- lm(prestige ~ income + education + type, data=Duncan)
+confidenceEllipses(m)
+confidenceEllipses(m, fill=TRUE)
+
+confidenceEllipses(Toy.mlm, fill=TRUE)
+
 
 
 #car::dfbetasPlots(Toy.lm1)

extra/confidenceEllipses.R

@@ -0,0 +1,95 @@
+confidenceEllipses <- function(model, ...) {
+  UseMethod("confidenceEllipses")
+}
+
+
+confidenceEllipses.default <-
+  function(model, coefnames,  main, ...) {
+    if (missing(main))
+      main <- paste("Pairwise Confidence Ellipses for",
+                    deparse(substitute(model)))
+    b <- coef(model)
+    p <- length(b)
+    if (missing(coefnames))
+      coefnames <- paste0(names(b), "\ncoefficient")
+    save <-
+      par(
+        mfrow = c(p, p),
+        mar = c(2, 2, 0, 0),
+        oma = c(0, 0, 2, 0)
+      )
+    on.exit(par(save))
+    ylab <- coefnames[1]
+    for (i in 1:p) {
+      for (j in 1:p) {
+        if (j == 1) {
+          yaxis <- TRUE
+        } else {
+          yaxis <- FALSE
+        }
+        if (i == p) {
+          xaxis <- TRUE
+        } else {
+          xaxis <- FALSE
+        }
+        if (i == j) {
+          if (i == 1) {
+            confidenceEllipse(
+              model,
+              c(2, 1),
+              xaxt = "n",
+              yaxt = "n",
+              center.pch = "",
+              col = "white",
+              grid = FALSE
+            )
+            axis(2)
+          } else if (j == p) {
+            confidenceEllipse(
+              model,
+              c(p, 2),
+              xaxt = "n",
+              yaxt = "n",
+              center.pch = "",
+              col = "white",
+              grid = FALSE
+            )
+            axis(1)
+          }
+          else {
+            confidenceEllipse(
+              model,
+              c(1, 2),
+              xaxt = "n",
+              yaxt = "n",
+              center.pch = "",
+              col = "white",
+              grid = FALSE
+            )
+          }
+          usr <- par("usr")
+          text(mean(usr[1:2]), mean(usr[3:4]), coefnames[i])
+        }
+        else{
+          confidenceEllipse(model, c(j, i), # xlab = xlab, ylab = ylab,
+                            xaxt = "n", yaxt = "n", ...)
+          if (j == 1)
+            axis(2)
+          if (i == p)
+            axis(1)
+        }
+      }
+    }
+    title(main = main,
+          outer = TRUE,
+          line = 1)
+    invisible(NULL)
+  }
+
+confidenceEllipses.mlm <- function(model, coefnames, ...) {
+  if (missing(coefnames))  {
+    coefnames <- rownames(vcov(model))
+    coefnames <- paste0(coefnames, "\ncoefficient")
+  }
+  NextMethod(coefnames = coefnames)
+}

extra/test-influencePlot.R

@@ -3,15 +3,17 @@ library(car)
 data(Duncan, package="carData")
 influencePlot(lm(prestige ~ income + education, data=Duncan))
 
-influencePlot(lm(prestige ~ income + education, data=Duncan), cex=2)
-
+# influencePlot(lm(prestige ~ income + education, data=Duncan), cex=2)
 # Error in plot.xy(xy.coords(x, y), type = type, ...) : 
 #   formal argument "cex" matched by multiple actual arguments
 
 # this works
-influencePlot(lm(prestige ~ income + education, data=Duncan), id=list(cex=2))
+influencePlot(lm(prestige ~ income + education, data=Duncan), id=list(cex=1.2))
 
 # test my modification to use fill ~ CookD
-applyDefaults <- car:::applyDefaults
-source("C:/R/Rprojects/mvinfluence/extra/influencePlot.R")
-influencePlot(lm(prestige ~ income + education, data=Duncan), id=list(cex=2))
+# applyDefaults <- car:::applyDefaults
+#source("C:/R/Rprojects/mvinfluence/extra/influencePlot.R")
+source(here::here("extra", "influencePlot.R"))
+influencePlot(lm(prestige ~ income + education, data=Duncan), id=list(cex=1.2))
+
+influencePlot(lm(prestige ~ income + education, data=Duncan), id=list(cex=1.2), fill.col = "red")