contrast.rms: made SE a vector not a matrix, added 4 list logic for n…

…vary, added new test from JoAnn Alvarez

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: rms
-Version: 4.3-1
-Date: 2015-04-20
+Version: 4.3-2
+Date: 2015-??-??
 Title: Regression Modeling Strategies
 Author: Frank E Harrell Jr <f.harrell@vanderbilt.edu>
 Maintainer: Frank E Harrell Jr <f.harrell@vanderbilt.edu>

NEWS

@@ -1,3 +1,6 @@
+Changes in version 4.3-2 (2015-??-??)
+   * contrast.rms: made SE a vector not a matrix, added 4 list logic for nvary, added new test from JoAnn Alvarez
+
 Changes in version 4.3-1 (2015-04-20)
    * NAMESPACE: removed reference to gridExtra, in DESCRIPTION moved gridExtra from Depends to Suggests
    * ggplot.Predict: re-worked generation of ggplot function call construction to use character strings with evaluation at the very end; added colfill argument

R/contrast.s

@@ -10,7 +10,6 @@ contrast.rms <-
   type <- match.arg(type)
   conf.type <- match.arg(conf.type)
   boot.type <- match.arg(boot.type)
-##  if(conf.type == 'simultaneous') require(multcomp)
 
   zcrit <- if(length(idf <- fit$df.residual)) qt((1 + conf.int) / 2, idf) else
               qnorm((1 + conf.int) / 2)
@@ -82,9 +81,11 @@ contrast.rms <-
   } else if(max(mall) > 1) {
     ## Label contrasts by values of longest variable in list if
     ## it has the same length as the expanded design matrix
-    d <- if(ma > 1) a else b  # TODO add logic for a2 b2
+    d <- if(ma > 1) a else b
+    if(! missing(a2) && (max(ma2, mb2) > max(ma, mb)))
+      d <- if(ma2 > 1) a2 else b2
     l <- sapply(d, length)
-    vary <- if(sum(l == max(ma, mb)) == 1) d[l == max(ma, mb)]
+    vary <- if(sum(l == max(mall)) == 1) d[l == max(mall)]
   }
 
   if(sum(mall > 1) > 1 && ! allsame(mall[mall > 1]))
@@ -117,8 +118,8 @@ contrast.rms <-
 
   est <- matxv(X, coef(fit))
   v <- X %*% vcov(fit, regcoef.only=FALSE) %*% t(X)
-  ndf <- if(is.matrix(v))nrow(v) else 1
-  se <- if(ndf==1) sqrt(v) else sqrt(diag(v))
+  ndf <- if(is.matrix(v)) nrow(v) else 1
+  se <- as.vector(if(ndf == 1) sqrt(v) else sqrt(diag(v)))
   Z <- est / se
   P <- if(length(idf)) 2 * (1 - pt(abs(Z), idf)) else 2 * (1 - pnorm(abs(Z)))
   if(conf.type != 'simultaneous') {
@@ -184,8 +185,9 @@ print.contrast.rms <- function(x, X=FALSE, fun=function(u) u,
   no[no=='Pvalue'] <- pn
 
   cnames <- x$cnames
-  if(! length(cnames)) cnames <- if(x$nvary) rep('', length(x[[1]])) else
-    as.character(1:length(x[[1]]))
+  if(! length(cnames))
+    cnames <- if(x$nvary) rep('', length(x[[1]])) else
+                          as.character(1 : length(x[[1]]))
   if(any(x$redundant)) cnames <- paste(ifelse(x$redundant, '*', ' '), cnames)
   w <- data.frame(w, row.names=paste(format(1:length(cnames)), cnames, sep=''))
   w$Contrast <- fun(w$Contrast)

tests/contrast.r

@@ -0,0 +1,33 @@
+# From JoAnn Alvarez
+
+require(rms)
+set.seed(1)
+age <- rnorm(200,40,12)
+sex <- factor(sample(c('female','male'),200,TRUE))
+logit <- (sex=='male') + (age-40)/5
+y <- ifelse(runif(200) <= plogis(logit), 1, 0)
+f <- lrm(y ~ pol(age,2)*sex)
+anova(f)
+# Compare a 30 year old female to a 40 year old male
+# (with or without age x sex interaction in the model)
+contrast(f, list(sex='female', age=30), list(sex='male', age=40))
+# Test for interaction between age and sex, duplicating anova
+contrast(f, list(sex='female', age=30),
+            list(sex='male',   age=30),
+            list(sex='female', age=c(40,50)),
+            list(sex='male',   age=c(40,50)), type='joint')
+# Duplicate overall sex effect in anova with 3 d.f.
+contrast(f, list(sex='female', age=c(30,40,50)),
+            list(sex='male',   age=c(30,40,50)), type='joint')
+
+
+jim <- contrast(f, list(sex = "male", age=30),
+                   list(sex = "male", age=40))
+print(jim, fun = exp)
+jane <- contrast(f, list(sex = c("male", "female"), age=30),
+                    list(sex = c("male", "female"), age=40))
+print(jane, fun = exp)
+
+
+
+