r - clmm后如何使用texreg(我想提取随机效应成分)

Question

在听取了@PhilipLeifeld 的建议(请参阅下面的评论部分)后，我正在根据自己的进步重写这篇博文。

我尝试使用 texreg 将 clmm 输出到 Latex。由于该包在其默认模式下不支持 clmm，因此我尝试使用 extract 功能扩展该包(参见 Print "beautiful" tables for h2o models in R 上的答案部分)。同时，我发现代码发布在 https://gist.github.com/kjgarza/340201f6564ca941fe25 上可以作为我的起点；我将代码称为下面的基线代码。以下模型(结果)非常代表我的实际代码。

library(ordinal)
library(texreg)
d<-data.frame(wine)
result<-clmm(rating~ 1+temp+contact+(1+temp|judge), data=d)

我想在 latex 表中显示的是随机效应成分，基线代码中省略了这些成分。以下是摘要输出的一部分。

summary(result)

Random effects:
 Groups Name        Variance Std.Dev. Corr  
 judge  (Intercept) 1.15608  1.0752         
        tempwarm    0.02801  0.1674   0.649 
 Number of groups:  judge 9 

具体来说，我想显示方差(和组数)；我不需要相关部分。在处理基线代码时，我还了解到“texreg”仅允许一组有限的 latex 显示参数，并且“include.variance”选项与我的目标相关。因此，我尝试将随机效应组件添加到“gof”参数，因为在基线代码中包含“include.variance”选项。

这是我所做的。首先，我在定义 extract.clmm 函数的部分添加了“include.variance”。

extract.clmm <- function(model, include.thresholds = TRUE, include.aic = TRUE, 
                     include.bic = TRUE, include.loglik = TRUE, include.variance = TRUE, oddsratios = TRUE, conf.level= 0.95, include.nobs = TRUE, ...) {
s <- summary(model, ...)

tab <- s$coefficients
thresh <- tab[rownames(tab) %in% names(s$alpha), ]
threshold.names <- rownames(thresh)
threshold.coef <- thresh[, 1]
threshold.se <- thresh[, 2]
threshold.pval <- thresh[, 4]
beta <- tab[rownames(tab) %in% names(s$beta), ]
beta.names <- rownames(beta)
beta.coef <- beta[, 1]
beta.se <- beta[, 2]
beta.pval <- beta[, 4]

然后，我添加了以下三行。

### for random effect components###   
rand<-s$ST[[1]]
rand.names<-rownames(rand)
rand.var<-rand[,1]

以下部分是我额外包含在基线代码中的部分(“include.variance”)。

if (include.variance == TRUE) {       
    gof.names <- c(gof.names, rand.names)
    gof <- c(gof, rand)
    gof.decimal <- c(gof.decimal, TRUE)   
}

运行 extract.clmm 函数后，我运行了以下命令。

test<-extract.clmm(result, include.variance=TRUE, oddsratios=FALSE)

然后，我收到一条错误消息:validityMethod(object) 错误:gof.names 和 gof 必须具有相同的长度! 而我发现“rand”和“rand. names”在“result”的情况下是 4 和 2，我不知道如何处理这个。任何意见将不胜感激。提前致谢。

Answer 1

让我们首先重写您的测试用例，使其同时包含一个具有随机效应的模型 (clmm) 和一个没有随机效应的模型 (clm)，两者均来自 ordinal。包裹。这将使我们能够检查 extract.clmm 是否我们将要编写的函数会产生以与现有 extract.clm 兼容的方式格式化的结果。 texreg 中的函数包裹:

library("ordinal")
library("texreg")
d <- data.frame(wine)
result.clmm <- clmm(rating ~ 1 + temp + contact + (1 + temp|judge), data = d)
result.clm <- clm(rating ~ 1 + temp + contact, data = d)

现有clm通用方法 extract在 texreg 中发挥作用看起来像这样，我们可以将其用作编写 clmm 的模板方法，因为两种对象类型的结构相似:

# extension for clm objects (ordinal package)
extract.clm <- function(model, include.thresholds = TRUE, include.aic = TRUE, 
    include.bic = TRUE, include.loglik = TRUE, include.nobs = TRUE, ...) {
  s <- summary(model, ...)

  tab <- s$coefficients
  thresh <- tab[rownames(tab) %in% names(s$aliased$alpha), , drop = FALSE]
  threshold.names <- rownames(thresh)
  threshold.coef <- thresh[, 1]
  threshold.se <- thresh[, 2]
  threshold.pval <- thresh[, 4]
  beta <- tab[rownames(tab) %in% names(s$aliased$beta), , drop = FALSE]
  beta.names <- rownames(beta)
  beta.coef <- beta[, 1]
  beta.se <- beta[, 2]
  beta.pval <- beta[, 4]
  if (include.thresholds == TRUE) {
    names <- c(beta.names, threshold.names)
    coef <- c(beta.coef, threshold.coef)
    se <- c(beta.se, threshold.se)
    pval <- c(beta.pval, threshold.pval)
  } else {
    names <- beta.names
    coef <- beta.coef
    se <- beta.se
    pval <- beta.pval
  }

  n <- nobs(model)
  lik <- logLik(model)[1]
  aic <- AIC(model)
  bic <- BIC(model)
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.aic == TRUE) {
    gof <- c(gof, aic)
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.bic == TRUE) {
    gof <- c(gof, bic)
    gof.names <- c(gof.names, "BIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.loglik == TRUE) {
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.nobs == TRUE) {
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num.\ obs.")
    gof.decimal <- c(gof.decimal, FALSE)
  }

  tr <- createTexreg(
      coef.names = names, 
      coef = coef, 
      se = se, 
      pvalues = pval, 
      gof.names = gof.names, 
      gof = gof, 
      gof.decimal = gof.decimal
  )
  return(tr)
}

setMethod("extract", signature = className("clm", "ordinal"), 
    definition = extract.clm)

clmm 的第一个差值objects 是系数​​等未存储在 summary(model)$aliased$alpha 下和 summary(model)$aliased$beta , 但直接在 summary(model)$alpha 之下和 summary(model)$beta .

我们需要做的第二件事是为组数和随机方差添加拟合优度元素。

组数显然存储在 summary(model)$dims$nlev.gf 下，具有针对不同条件变量的多个条目。所以这很容易。

随机方差没有存储在任何地方，所以我们需要在 source code of the ordinal package 中查找它.我们可以在那里看到 print.summary.clmm函数使用名为 formatVC 的内部辅助函数打印差异。此函数包含在相同的 R 中脚本，基本上只是进行格式化并调用另一个名为 varcov 的内部辅助函数(也包含在同一个文件中)来计算方差。该函数反过来计算 model$ST 的转置叉积得到方差。我们可以直接在 extract.clmm 的 GOF block 中简单地做同样的事情。功能，例如，使用 diag(s$ST[[1]] %*% t(s$ST[[1]]))对于第一个随机效应。我们只需要确保对所有随机效果都这样做，这意味着我们需要将其放入循环中并替换 [[1]]通过像 [[i]] 这样的迭代器.

决赛clmm extract 的方法函数可能如下所示:

# extension for clmm objects (ordinal package)
extract.clmm <- function(model, include.thresholds = TRUE,
    include.loglik = TRUE, include.aic = TRUE,  include.bic = TRUE,
    include.nobs = TRUE, include.groups = TRUE, include.variance = TRUE, ...) {
  s <- summary(model, ...)

  tab <- s$coefficients
  thresh <- tab[rownames(tab) %in% names(s$alpha), ]
  threshold.names <- rownames(thresh)
  threshold.coef <- thresh[, 1]
  threshold.se <- thresh[, 2]
  threshold.pval <- thresh[, 4]
  beta <- tab[rownames(tab) %in% names(s$beta), ]
  beta.names <- rownames(beta)
  beta.coef <- beta[, 1]
  beta.se <- beta[, 2]
  beta.pval <- beta[, 4]

  if (include.thresholds == TRUE) {
    cfnames <- c(beta.names, threshold.names)
    coef <- c(beta.coef, threshold.coef)
    se <- c(beta.se, threshold.se)
    pval <- c(beta.pval, threshold.pval)
  } else {
    cfnames <- beta.names
    coef <- beta.coef
    se <- beta.se
    pval <- beta.pval
  }

  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.loglik == TRUE) {
    lik <- logLik(model)[1]
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.aic == TRUE) {
    aic <- AIC(model)
    gof <- c(gof, aic)
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.bic == TRUE) {
    bic <- BIC(model)
    gof <- c(gof, bic)
    gof.names <- c(gof.names, "BIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.nobs == TRUE) {
    n <- nobs(model)
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num.\ obs.")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  if (include.groups == TRUE) {
    grp <- s$dims$nlev.gf
    grp.names <- paste0("Groups (", names(grp), ")")
    gof <- c(gof, grp)
    gof.names <- c(gof.names, grp.names)
    gof.decimal <- c(gof.decimal, rep(FALSE, length(grp)))
  }
  if (include.variance == TRUE) {
    var.names <- character()
    var.values <- numeric()
    for (i in 1:length(s$ST)) {
      variances <- diag(s$ST[[i]] %*% t(s$ST[[i]]))
      var.names <- c(var.names, paste0("Variance: ", names(s$ST)[[i]], ": ", 
          names(variances)))
      var.values <- c(var.values, variances)
    }
    gof <- c(gof, var.values)
    gof.names <- c(gof.names, var.names)
    gof.decimal <- c(gof.decimal, rep(TRUE, length(var.values)))
  }

  tr <- createTexreg(
      coef.names = cfnames, 
      coef = coef, 
      se = se, 
      pvalues = pval, 
      gof.names = gof.names, 
      gof = gof, 
      gof.decimal = gof.decimal
  )
  return(tr)
}

setMethod("extract", signature = className("clmm", "ordinal"), 
    definition = extract.clmm)

您可以只在运行时执行代码，texreg应该能够从 clmm 创建表对象，包括随机方差。我将把这段代码添加到下一个 texreg发布。

您可以将其应用到您的示例中，如下所示:

screenreg(list(result.clmm, result.clm), single.row = TRUE)

结果在 clmm 之间是兼容的和 clm对象，正如您在输出中看到的那样:

==================================================================
                              Model 1            Model 2          
------------------------------------------------------------------
tempwarm                        3.07 (0.61) ***    2.50 (0.53) ***
contactyes                      1.83 (0.52) ***    1.53 (0.48) ** 
1|2                            -1.60 (0.69) *     -1.34 (0.52) ** 
2|3                             1.50 (0.60) *      1.25 (0.44) ** 
3|4                             4.22 (0.82) ***    3.47 (0.60) ***
4|5                             6.11 (1.02) ***    5.01 (0.73) ***
------------------------------------------------------------------
Log Likelihood                -81.55             -86.49           
AIC                           181.09             184.98           
BIC                           201.58             198.64           
Num. obs.                      72                 72              
Groups (judge)                  9                                 
Variance: judge: (Intercept)    1.16                              
Variance: judge: tempwarm       0.03                              
==================================================================
*** p < 0.001, ** p < 0.01, * p < 0.05

您可以使用参数 include.variances == FALSE和 include.groups == FALSE如果需要，可以关闭方差和组大小的报告。