harrelfe / rms

Regression Modeling Strategies
https://hbiostat.org/R/rms
Other
170 stars 48 forks source link

print.Gls() using corSymm() returns error message. #136

Open SatsukiTaniuchi opened 9 months ago

SatsukiTaniuchi commented 9 months ago

Dear Prof. Harrell,

Gls() using corSymm() in rms package version 6.8-0 returns an error message when displaying the results of the fit. The results of the fit can be displayed when using compound symmetry correlation.

library(nlme)
library(rms)        # remotes::install_github("harrelfe/rms")
library(dplyr)
library(tidyverse)
# _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
# create data
# _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
beta0   = 100
betag   = 10
beta1   = 20
beta2   = 30
betag1  = 40
betag2  = 50
vmat    = matrix(c(400, 200, 150, 200, 500, 100, 150, 100, 400), nrow = 3)
mu0     = c(beta0        , beta0 + beta1                 , beta0 + beta2                 )
mu1     = c(beta0 + betag, beta0 + beta1 + betag + betag1, beta0 + beta2 + betag + betag2)

set.seed(1)

data = rbind( 
     data.frame(id =   1:100, group = "C", MASS::mvrnorm(n = 100, mu = mu0, Sigma = vmat))
    ,data.frame(id = 101:200, group = "T", MASS::mvrnorm(n = 100, mu = mu1, Sigma = vmat))
)
colnames(data) = c("id", "group", "y0", "y1", "y2")
data = data %>%
       pivot_longer(cols = matches("y\\d"), names_pattern = "(y)(\\d)", names_to = c(".value","time")) %>%
       mutate(across(
            .cols = c(time, group)
           ,.fns  = ~ factor(.x)  
       ))
# _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
# fitting
# _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
fit.nlme = gls( model       = y ~ group * time
               ,data        = data
               ,correlation = nlme::corSymm(form = ~ as.numeric(time)|id)
               ,weights     = nlme::varIdent(form = ~ 1 | time)
               ,na.action   = na.exclude
               ,method      = "REML"
           )
fit.nlme
print(fit.nlme) # displayed

fit.rms = Gls( model       = y ~ group * time
              ,data        = data
              ,correlation = nlme::corSymm(form = ~ as.numeric(time)|id)
              ,weights     = nlme::varIdent(form = ~ 1 | time)
              ,na.action   = na.exclude
              ,method      = "REML"
          )

# old version(rms:6.0-0, Hmisc:4.4-0) -> I can print result of fitting. 
# latest version(rms:6.8-0, Hmisc:5.1-1) -> Fitting  successed but I can't print result of fitting. 
contrast(fit.rms, list(time = "1", group = "T"), list(time = "1", group = "C"))

fit.rms
print(fit.rms)
#Error in print.default(x[-1, -p, drop = FALSE], ..., quote = FALSE) : 
#  formal argument "quote" matched by multiple actual arguments

traceback()
#14: print(x[-1, -p, drop = FALSE], ..., quote = FALSE)
#13: print.correlation(val, ...)
#12: print(val, ...)
#11: print.summary.corSymm(X[[i]], ...)
#10: FUN(X[[i]], ...)
#9: lapply(x, print, ...)

#... omitted

#6: do.call(type, c(obj, list(quote = FALSE)))
#5: withVisible(...elt(i))
#4: capture.output(do.call(type, c(obj, list(quote = FALSE))))
#3: prModFit(x, title = title, z, digits = digits, coefs = coefs, 
#       ...)
#2: print.Gls(fit.rms)
#1: print(fit.rms)

fit.rms = Gls( model       = y ~ group * time
              ,data        = data
              ,correlation = nlme::corCompSymm(form = ~ as.numeric(time)|id)
              ,na.action   = na.exclude
              ,method      = "REML"
          )
print(fit.rms) # displayed

My apologies for the redundant example.