clmm后如何使用texreg（我想提取随机效应成分）
How to use texreg after clmm (I want to extract random effect components)
我正在根据@PhilipLeifeld 的建议（请参阅下面的评论部分）根据我的进度重写此帖子。
我尝试使用 texreg 将 clmm 输出到 Latex。由于包在默认模式下不支持 clmm，我尝试使用 extract 功能扩展包（请参阅 ). Meanwhile, I found that code posted on 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 table 中显示的是 运行dom 效果组件，这些组件在基线代码中被省略。以下是摘要输出的一部分。
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" 只允许一组有限的乳胶显示参数，并且 "include.variance" 的选项与我的目标相关。因此，我尝试将 运行dom 效果组件添加到 "gof" 参数，因为在基线代码中包含 "include.variance" 选项。
这是我所做的。首先，我将 "include.variance" 添加到定义 extract.clmm 函数的部分。
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，我不知道如何处理。任何意见将不胜感激。提前致谢。
让我们首先重写您的测试用例，使其同时包含一个具有随机效应的模型 (clmm) 和一个没有随机效应的模型 (clm)，两者均来自 ordinal包裹。这将允许我们检查我们将要编写的 extract.clmm 函数生成的结果是否以与 texreg 包中现有的 extract.clm 函数兼容的方式格式化：
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)
texreg 中通用 extract 函数的现有 clm 方法如下所示，我们可以将其用作编写 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 对象的第一个区别是系数等不存储在 summary(model)$aliased$alpha 和 summary(model)$aliased$beta 下，而是直接存储在 summary(model)$alpha 和 [=29= 下].
我们需要做的第二件事是为组数和随机方差添加拟合优度元素。
组数显然存储在 summary(model)$dims$nlev.gf 下，其中包含针对不同条件变量的多个条目。所以很简单。
随机方差没有存储在任何地方，所以我们需要在 source code of the ordinal package 中查找它。我们可以看到 print.summary.clmm 函数使用一个名为 formatVC 的内部辅助函数来打印方差。此函数包含在同一个 R 脚本中，基本上只是进行格式化并调用另一个名为 varcov 的内部辅助函数（也包含在同一文件中）来计算方差。该函数反过来计算 model$ST 的转置叉积以获得方差。我们可以直接在 extract.clmm 函数的 GOF 块中直接做同样的事情，例如，使用 diag(s$ST[[1]] %*% t(s$ST[[1]])) 作为第一个随机效应。我们只需要确保对所有随机效果都这样做，这意味着我们需要将其放入循环中并用 [[i]].
这样的迭代器替换 [[1]]
extract 函数的最终 clmm 方法可能如下所示：
# 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 关闭方差和组大小的报告。
作为对@Philip 的回答的快速评论，在新版本或 R studio 中，以下内容不 return 矩阵：
thresh <- tab[rownames(tab) %in% names(s$alpha), ]
这会导致以下代码 return 出错。但是，对此的快速修复可以是：
thresh <- subset.matrix(tab, rownames(tab) %in% names(s$alpha) )
beta <- subset.matrix(tab, rownames(tab) %in% names(s$beta) )
我正在根据@PhilipLeifeld 的建议（请参阅下面的评论部分）根据我的进度重写此帖子。
我尝试使用 texreg 将 clmm 输出到 Latex。由于包在默认模式下不支持 clmm，我尝试使用 extract 功能扩展包（请参阅
library(ordinal)
library(texreg)
d<-data.frame(wine)
result<-clmm(rating~ 1+temp+contact+(1+temp|judge), data=d)
我想在 Latex table 中显示的是 运行dom 效果组件，这些组件在基线代码中被省略。以下是摘要输出的一部分。
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" 只允许一组有限的乳胶显示参数，并且 "include.variance" 的选项与我的目标相关。因此，我尝试将 运行dom 效果组件添加到 "gof" 参数，因为在基线代码中包含 "include.variance" 选项。
这是我所做的。首先，我将 "include.variance" 添加到定义 extract.clmm 函数的部分。
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，我不知道如何处理。任何意见将不胜感激。提前致谢。
让我们首先重写您的测试用例，使其同时包含一个具有随机效应的模型 (clmm) 和一个没有随机效应的模型 (clm)，两者均来自 ordinal包裹。这将允许我们检查我们将要编写的 extract.clmm 函数生成的结果是否以与 texreg 包中现有的 extract.clm 函数兼容的方式格式化：
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)
texreg 中通用 extract 函数的现有 clm 方法如下所示，我们可以将其用作编写 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 对象的第一个区别是系数等不存储在 summary(model)$aliased$alpha 和 summary(model)$aliased$beta 下，而是直接存储在 summary(model)$alpha 和 [=29= 下].
我们需要做的第二件事是为组数和随机方差添加拟合优度元素。
组数显然存储在 summary(model)$dims$nlev.gf 下，其中包含针对不同条件变量的多个条目。所以很简单。
随机方差没有存储在任何地方，所以我们需要在 source code of the ordinal package 中查找它。我们可以看到 print.summary.clmm 函数使用一个名为 formatVC 的内部辅助函数来打印方差。此函数包含在同一个 R 脚本中，基本上只是进行格式化并调用另一个名为 varcov 的内部辅助函数（也包含在同一文件中）来计算方差。该函数反过来计算 model$ST 的转置叉积以获得方差。我们可以直接在 extract.clmm 函数的 GOF 块中直接做同样的事情，例如，使用 diag(s$ST[[1]] %*% t(s$ST[[1]])) 作为第一个随机效应。我们只需要确保对所有随机效果都这样做，这意味着我们需要将其放入循环中并用 [[i]].
[[1]]
extract 函数的最终 clmm 方法可能如下所示：
# 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 关闭方差和组大小的报告。
作为对@Philip 的回答的快速评论，在新版本或 R studio 中，以下内容不 return 矩阵：
thresh <- tab[rownames(tab) %in% names(s$alpha), ]
这会导致以下代码 return 出错。但是，对此的快速修复可以是：
thresh <- subset.matrix(tab, rownames(tab) %in% names(s$alpha) )
beta <- subset.matrix(tab, rownames(tab) %in% names(s$beta) )