tests/effect-tests-2.R

In effects: Effect Displays for Linear, Generalized Linear, and Other Models

if (requireNamespace("carData") && require("effects")){

  # plots should show fitted values directly on plotted effect, and must be checked visually
  # numbering corresponds to effect-test-1.R
  
  data(Duncan, package="carData")
  
  mod.1 <- lm(prestige ~ type + poly(income, degree=2, raw=TRUE), data=Duncan)
  

  # (2) focal: polynomial, constant: factor
  
  print(plot(Effect(c("income"), mod.1, residuals=TRUE), show.fitted=TRUE))
  
  if (!isTRUE(all.equal(Effect(c("income"), mod.1, residual=TRUE)$fit, 
                        Effect(c("income"), mod.1, 
                               xlevels=list(income=seq(7, 81, length.out=100)))$fit)))
    stop("failed test 2 (2)")
  
  
  # (3) focal: factor*polynomial, constant: polynomial
  
  mod.2 <- lm(prestige ~ type*poly(income, degree=2, raw=TRUE) + 
                poly(education, degree=2, raw=TRUE), data=Duncan)

  print(plot(Effect(c("type", "income"), mod.2, residuals=TRUE), show.fitted=TRUE))
  
  if (!isTRUE(all.equal(Effect(c("type", "income"), mod.2, residuals=TRUE)$fit, 
                        Effect(c("type", "income"), mod.2, 
                               xlevels=list(income=seq(7, 81, length.out=100)))$fit)))
    stop("failed test 2 (3)")
  
  
  # (4) focal: polynomial, constant: factor*polynomial
  
  print(plot(Effect(c("education"), mod.2, residuals=TRUE), show.fitted=TRUE))
  
  if (!isTRUE(all.equal(Effect(c("education"), mod.2, residuals=TRUE)$fit, 
                        Effect(c("education"), mod.2, 
                               xlevels=list(education=seq(7, 100, length.out=100)))$fit)))
    stop("failed test 2 (4)")  
  
  
  # (6) focal: factor*polynomial, constant: polynomial, using predict() & orthog. polys.
  
  mod.3 <- lm(prestige ~ type*poly(income, degree=2) + poly(education, degree=2), data=Duncan)
  
  print(plot(Effect(c("type", "income"), mod.3, residuals=TRUE), show.fitted=TRUE))
  
  if (!isTRUE(all.equal(Effect(c("type", "income"), mod.3, residuals=TRUE)$fit, 
                        Effect(c("type", "income"), mod.3, 
                               xlevels=list(income=seq(7, 81, length.out=100)))$fit)))
    stop("failed test 2 (6)")  
  
  # (7) focal: factor, constant: poly*poly
  
  mod.4 <- lm(prestige ~ type + poly(income, 2)*poly(education, 2), data=Duncan)
  
  print(plot(Effect(c("income", "education"), mod.4, residuals=TRUE), show.fitted=TRUE))
  
  if (!isTRUE(all.equal(Effect(c("income", "education"), mod.4, residuals=TRUE)$fit, 
                        Effect(c("income", "education"), mod.4, 
                               xlevels=list(income=seq(7, 81, length.out=100), 
                                            education=quantile(Duncan$education, probs=seq(0.2, 0.8, by=0.2))))$fit)))
    stop("failed test 2 (7)")  
  
  # (9) focal: covariate, constant: 2 factors and 1 covariate, 3-way interaction
  
  data(Mroz, package="carData")
  mod.6 <- lm(lwg ~ inc + age*hc*wc, data=Mroz)
  inc <- range(Mroz$inc)
  age <- range(Mroz$age)

  print(plot(Effect(c("inc"), mod.6, residuals=TRUE), show.fitted=TRUE))
  
  if (!isTRUE(all.equal(Effect(c("inc"), mod.6, residuals=TRUE)$fit, 
                        Effect(c("inc"), mod.6, 
                               xlevels=list(inc=seq(inc[1], inc[2], length.out=100)))$fit)))
    stop("failed test 2 (9-1)")  
  
  print(plot(Effect(c("age", "hc", "wc"), mod.6, residuals=TRUE), show.fitted=TRUE))
  
  if (!isTRUE(all.equal(Effect(c("age", "hc", "wc"), mod.6, residuals=TRUE)$fit, 
                        Effect(c("age", "hc", "wc"), mod.6, 
                               xlevels=list(age=seq(age[1], age[2], length.out=100)))$fit)))
    stop("failed test 2 (9-2)")  
  
  # additional tests of partial residuals
  
  income <- range(na.omit(Prestige)$income)
  
  mod.7 <- lm(prestige ~ income*type + education, data=Prestige)
  print(plot(Effect(c("income", "type"), mod.7, residuals=TRUE), show.fitted=TRUE))
  
  if (!isTRUE(all.equal(Effect(c("income", "type"), mod.7, residuals=TRUE)$fit, 
                        Effect(c("income", "type"), mod.7, 
                               xlevels=list(income=seq(income[1], income[2], length.out=100)))$fit)))
    stop("failed test 2 (additional-1)")  
  
  Mroz2 <- Mroz
  Mroz2$hc <- as.numeric(Mroz$hc) - 1
  Mroz2$wc <- as.numeric(Mroz$wc) - 1
  inc <- range(Mroz2$inc)
  mod.8 <- lm(lwg ~ inc*age*k5 + hc*wc, data=Mroz2)
  print(plot(Effect(c("inc", "age", "k5"), mod.8, residuals=TRUE, xlevels=list(k5=0:1)), 
             show.fitted=TRUE))
  
  if (!isTRUE(all.equal(Effect(c("inc", "age", "k5"), mod.8, residuals=TRUE, xlevels=list(k5=0:1))$fit, 
                        Effect(c("inc", "age", "k5"), 
                               mod.8, residuals=TRUE, 
                               xlevels=list(k5=0:1, inc=seq(inc[1], inc[2], length.out=100), age=quantile(Mroz2$age, seq(.2, .8, by=.2))))$fit)))
    stop("failed test 2 (additional-2)")  
  
  print(plot(Effect(c("hc", "wc"), mod.8, residuals=TRUE, xlevels=list(hc=0:1, wc=0:1)), 
       show.fitted=TRUE, smooth.residuals=FALSE, residuals.pch="."))
}