Убедитесь, что модель имеет только один фактор-ковариат - PullRequest
Новая
       83

Убедитесь, что модель имеет только один фактор-ковариат

0 голосов
/ 24 мая 2018

Я пишу пакет R, где основная функция принимает модель , которая может иметь только однофакторный ковариат (смещения допускаются).Чтобы убедиться, что пользователь соблюдает это правило, мне нужно проверить это.

В качестве примера рассмотрим следующие четыре модели: 

set.seed(123)
n <- 10 

## data
data <- data.frame(y = rnorm(n), 
  trt = rep(c(0, 1), each = n/2),
  x = 1:n)
datan <- data
datan$trt <- as.factor(datan$trt)

## models
mod1 <- lm(y ~ factor(trt), data = data)
mod2 <- lm(y ~ offset(x) + as.factor(trt), data = data)
mod3 <- lm(y ~ trt, data = datan)
mod4 <- glm(y ~ trt + offset(x), data = data)
mod5 <- lm(y ~ x + as.factor(trt), data = data)

Модели 1, 2 и 3 в порядке , модели 4 и 5 -не в порядке (модель 4 имеет нефакторную переменную trt, модель 5 имеет второй ковариат x).

Как проверить это с помощью R?Оптимально я бы получил TRUE для модели, которая в порядке, и FALSE для проблемной модели.

Это должно работать не только с lm() и glm(), но также с survreg() и coxph() (из пакета выживания).Кое-что, что может быть полезно, смотрит на формулу eval(getCall(mod1)$formula) и данные (data / datan).

Ответы [ 2 ]

0 голосов
/ 24 мая 2018

Как указано в предыдущем ответе @LAP, вы можете использовать terms() из этих моделей.Однако я бы порекомендовал взглянуть на attr(..., "factors") и attr(..., "dataClasses") вместо того, чтобы перейти к $model, который требует, чтобы в модели хранилось все model.frame().Это может или не может иметь место.В частности, при повторной подгонке нескольких моделей вы можете захотеть не сохранять фрейм модели каждый раз.

Таким образом, одной из идей будет продолжение следующих шагов:

  • Проверьте, не содержит ли attr(..., "factors") точно один столбец, и вы можете вернуть FALSE.
  • Если указан ровно один фактор, вы можете проверить соответствующий attr(..., "dataClasses"), если он равен "factor" / "ordered"и затем верните TRUE, в противном случае FALSE.

R код: 

one_factor <- function(object) {
  f <- attr(terms(object), "factors")
  if(length(f) == 0L || NCOL(f) != 1L) return(FALSE)
  d <- attr(terms(object), "dataClasses")
  if(d[colnames(f)] %in% c("ordered", "factor")) {
    return(TRUE)
  } else {
    return(FALSE)
  }
}

Похоже, это хорошо работает для объектов, основанных на одной части formula.

Фиктивные данные с числовым / коэффициентом / упорядоченным trt: 

d1 <- d2 <- d3 <- data.frame(y = log(1:9), x = 1:9, trt = rep(1:3, each = 3)) 
d2$trt <- factor(d2$trt)
d3$trt <- ordered(d3$trt)

Различные спецификации формул: 

f <- list(
  y ~ 1,
  y ~ x,
  y ~ trt,
  y ~ trt + x,
  y ~ trt + offset(x),
  y ~ trt + x + offset(x),
  y ~ trt + offset(as.numeric(trt)),
  y ~ factor(trt),
  y ~ factor(trt) + offset(x),
  y ~ factor(x > as.numeric(trt)),
  y ~ interaction(x, trt),
  y ~ 0 + trt
)

Ожидаемые результаты для d1, d2,и d3, соответственно: 

ok1 <- c(FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE)
ok2 <- c(FALSE, FALSE, TRUE,  FALSE, TRUE,  FALSE, TRUE,  TRUE, TRUE, TRUE, TRUE, TRUE)
ok3 <- ok2

Проверки для lm без сохранения кадра модели: 

lm1 <- lapply(f, lm, data = d1, model = FALSE)
identical(sapply(lm1, one_factor), ok1)
## [1] TRUE
lm2 <- lapply(f, lm, data = d2, model = FALSE)
identical(sapply(lm2, one_factor), ok2)
## [1] TRUE
lm3 <- lapply(f, lm, data = d3, model = FALSE)
identical(sapply(lm3, one_factor), ok3)
## [1] TRUE

Проверки для survreg (гауссовский) и coxph.(Последний выдает много предупреждений о несовпадении, что неудивительно, учитывая фиктивную структуру данных. Проверки по-прежнему работают, как и предполагалось.) 

library("survival")
d1$y <- d2$y <- d3$y <- Surv(d1$y + 0.5)

sr1 <- lapply(f, survreg, data = d1)
identical(sapply(sr1, one_factor), ok1)
## [1] TRUE
sr2 <- lapply(f, survreg, data = d2)
identical(sapply(sr2, one_factor), ok2)
## [1] TRUE
sr3 <- lapply(f, survreg, data = d3)
identical(sapply(sr3, one_factor), ok3)
## [1] TRUE

cph1 <- lapply(f, coxph, data = d1)
identical(sapply(cph1, one_factor), ok1)
## [1] TRUE
cph2 <- lapply(f, coxph, data = d2)
identical(sapply(cph2, one_factor), ok2)
## [1] TRUE
cph3 <- lapply(f, coxph, data = d3)
identical(sapply(cph3, one_factor), ok3)
## [1] TRUE

Примечание: Если у вас естьВ составных объектах Formula на основе этой функции может произойти сбой, и базовые тесты необходимо будет адаптировать.Примерами последних могут быть модели регрессии счета (zeroinfl, hurdle), полиномиальный логит (mlogit), инструментальные переменные (ivreg), гетероскедастические модели (vglm, betareg, crch)и т. д. У них могут быть такие формулы, как y ~ trt | 1 или y ~ trt | trt или y ~ trt | x, которые могут или не могут быть осуществимы в вашей структуре.

0 голосов
/ 24 мая 2018

Это потребует дополнительного тестирования, но оно работает для ваших примеров: 

FOO <- function(x){
  vars <- labels(terms(x))
  test <- sapply(x$model[vars], class)
  all(test == "factor", length(test) == 1)
}

Сначала мы извлекаем ковариаты модели, используя labels(terms()), что дает дополнительное преимущество, чтобы игнорировать смещения, затем получаемвектор классов и проверьте, выполняются ли два условия (1. переменная является фактором, 2. это только одна переменная). 

> sapply(list(mod1, mod2, mod3, mod4, mod5), FOO)
[1]  TRUE  TRUE  TRUE FALSE FALSE
...