Как использовать texreg после clmm (я хочу извлечь компоненты случайного эффекта)
Я переписываю эту публикацию, основываясь на моем прогрессе после получения советов от @PhilipLeifeld (см. Раздел комментариев ниже).
Я пытался поставить clmm выводит в латекс, используя texreg, Поскольку пакет не поддерживает clmm в режиме по умолчанию я попытался расширить пакет extract функция (см. часть ответа на Распечатка "красивых" таблиц для моделей H2O в 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)
То, что я хотел бы отобразить в латексной таблице, это компоненты случайных эффектов, которые опущены в базовом коде. Следующее является частью итогового результата.
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" имеет отношение к моей цели. Таким образом, я попытался добавить компоненты случайных эффектов в аргумент "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(объект): gof.names и gof должны иметь одинаковую длину! Хотя я обнаружил, что длины "rand" и "rand.names" в случае "result" равны 4 и 2, я не знаю, как с этим справиться. Любые комментарии будут очень благодарны. Заранее спасибо.
2 ответа
Давайте сначала перепишем ваш тестовый пример так, чтобы он содержал как модель со случайными эффектами (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 объекты в том, что коэффициенты и т. д. не хранятся под summary(model)$aliased$alpha а также summary(model)$aliased$beta, но прямо под summary(model)$alpha а также summary(model)$beta,
Второе, что нам нужно сделать, это добавить элементы соответствия для количества групп и случайных отклонений.
Количество групп, по-видимому, хранится в summary(model)$dims$nlev.gf, с несколькими записями для различных переменных кондиционирования. Так легко.
Случайные отклонения нигде не хранятся, поэтому мы должны посмотреть это в исходном коде ordinal пакет Там видно, что print.summary.clmm функция использует внутреннюю вспомогательную функцию под названием formatVC распечатать дисперсии. Эта функция содержится в том же R сценарий и в основном просто делает форматирование и вызывает другую внутреннюю вспомогательную функцию под названием varcov (также содержится в том же файле) для расчета отклонений. Эта функция, в свою очередь, вычисляет транспонированный перекрестный продукт model$ST чтобы получить отклонения. Мы можем просто сделать то же самое прямо в блоке GOF нашего extract.clmm функция, например, используя 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 чтобы отключить отчеты о дисперсиях и размерах групп, если хотите.
В качестве быстрого замечания по ответу @ Филиппа, в новой версии или R studio следующее не возвращает матрицу:
thresh <- tab[rownames(tab) %in% names(s$alpha), ]
Это приводит к тому, что следующий код возвращает ошибку. Однако, быстрое решение для этого может быть:
thresh <- subset.matrix(tab, rownames(tab) %in% names(s$alpha) )
beta <- subset.matrix(tab, rownames(tab) %in% names(s$beta) )