Использование функции прогнозирования на биномиальном GLMM: получение числового, а не биномиального значения - PullRequest
Новая
       58

Использование функции прогнозирования на биномиальном GLMM: получение числового, а не биномиального значения

0 голосов
/ 23 мая 2019

Я использовал GLMM, используя пакет glmmTMB для большого набора данных.Пример данных ниже 

structure(list(code = structure(c(1L, 1L, 1L, 1L), .Label = c("2388", 
"2390", "12950", "12952", "12954", "12956", "12958", "12960", 
"12962", "12964", "12966", "12968", "13573", "13574", "13575", 
"13576", "13577", "14203", "19318", "19319", "19320", "19321", 
"19322", "19515", "19517", "19523", "19524", "25534", "25535", 
"25536", "25537", "25539", "25540", "25541", "25542", "25543", 
"25544", "25545", "25546", "25547", "25548", "25549", "25550", 
"25552", "25553", "27583", "27584", "27585", "27586", "27588", 
"27589", "27590", "27591", "27592", "27593", "27594", "27595", 
"27596", "27597", "27598", "27599", "27600", "27601", "27602", 
"27603", "27604", "27605", "27606", "27607", "27608", "27609", 
"27610", "27611", "27613", "27614", "27615", "27616", "27617", 
"27618", "27619", "27620", "27621", "27622", "27624", "27625", 
"27626", "27627", "27629", "27630", "27631", "27632", "34176", 
"34177", "34178", "34179", "52975", "52977", "52978", "54814", 
"54815", "54816", "54817", "54821", "54822", "54823", "54824", 
"54825", "54835", "54837", "54838", "54839", "54840", "54841", 
"54842", "54843", "54844", "54845", "54846", "54847", "54848", 
"54849", "54851", "54852", "54853", "54856", "54858", "54859", 
"54860", "54863", "54864", "54866", "54869", "54872", "54873", 
"54874", "54875", "54876", "54877", "54878", "54880", "54882", 
"54883", "54884", "54886", "54887", "54889", "54890", "54892", 
"54893", "54895", "54896", "54898", "54899", "54900", "54901", 
"54903", "54904", "54905", "54906", "54911", "54912", "54914", 
"54915", "54931", "54933", "54934", "54935", "54937", "54939", 
"54940", "54941", "54942", "54943", "54944", "54945", "54946", 
"54947", "54948", "54950", "54952", "54954", "54955", "54957", 
"54958", "54959", "54961", "54962"), class = "factor"), station = 
c("PB14", 
"PB14", "PB16", "PB16"), species = c("Silvertip Shark", "Silvertip Shark", 
"Silvertip Shark", "Silvertip Shark"), sex = c("F", "F", "F", 
"F"), size = c(112, 112, 112, 112), datetime = c("1466247120", 
"1466247420", "1467026100", "1469621400"), year = c("2016", "2016", 
"2016", "2016"), month = c(6, 6, 6, 7), hour = c("11", "11", 
"12", "13"), season = c("dry season", "dry season", "dry season", 
"dry season"), daynight = c("day", "day", "day", "day"), time_diff = c(4, 
5, 5821, 43255), offshore = structure(c(2L, 2L, 1L, 1L), .Label = 
c("offshore", 
"onshore"), class = "factor"), rowN = 1:4), row.names = c(NA, 
4L), class = "data.frame")

Я собираюсь запустить модель на 80% моих данных, а затем проверить ее с помощью функции предсказания ().Используя приведенный ниже код 

Off_80 <- Off %>% sample_frac(.80)

Off_20  <- anti_join(Off, Off_80, by = 'rowN')

OffMod_80 <- glmmTMB(offshore ~ sex + log(size) + species*daynight + species*season + (1|code), family=binomial(), data=Off_80)

pred_Off20 <- as.data.frame(predict(OffMod_80, newdata=Off_20, type="response", allow.new.levels=TRUE))

Идея заключается в том, что я бы сравнил прогнозируемые результаты с моими наблюдаемыми результатами, чтобы проверить силу модели.'или' береговый 'ответ, я получаю числовое значение.

         predict()
     1   0.2807461                                                                     
     2   0.2631816                                                                       
     3   0.2631816                                                                     
     4   0.2807461

Есть ли способ заставить функцию предсказания выплевывать биномиальный ответ?Или как мне интерпретировать эти значения как биномиальные?

Первоначально у меня была переменная ответа как 1 или 0, но после этого post я изменил значения на факторы.Но по-прежнему предикат () показывает числовое значение.

Любая помощь приветствуется!

Ответы [ 2 ]

1 голос
/ 29 мая 2019

predict() для биномиальных моделей возвращает вероятность успеха или неудачи, но не возвращает 1 или 0 (поскольку вы можете предсказать этот результат только с определенной вероятностью ),Поэтому, если вы хотите проверить производительность вашей модели, вы можете попробовать вычислить область под кривой : 

library(glmmTMB)
library(pROC)

data(mtcars)
n <- nrow(mtcars)
train <- sample(1:n, n * .8, TRUE)
mtcars_train <- mtcars[train, ]
mtcars_test <- mtcars[-train, ]

m <- glmmTMB(formula = vs ~ hp + wt + (1 | gear), family = binomial, data = mtcars_train)
p <- predict(m, newdata = mtcars_test)

auc(roc(response = mtcars_test$vs, predictor = p))
#> Area under the curve: 0.9107

Создано в 2019-05-29 с помощью Представить пакет (v0.3.0)

Другой способ - сравнить процент правильных прогнозов (PCP) полной модели с нулевой моделью (см. Herron, M. (1999). Постоценка неопределенности в моделях с ограниченными зависимыми переменными. Политический анализ, 8, 83–98 ).Здесь PCP полной модели должен быть значительно выше: 

library(glmmTMB)
library(insight)

data(mtcars)

m <- glmmTMB(formula = vs ~ hp + wt + (1 | gear), family = binomial, data = mtcars)
m0 <- glmmTMB(formula = vs ~ 1 + (1 | gear), family = binomial, data = mtcars)

y <- insight::get_response(m)
y0 <- insight::get_response(m0)

n <- nobs(m)
n0 <- nobs(m0)

p <- predict(m, type = "response")
p0 <- predict(m0, type = "response")

pcp_full <- (sum(1 - p[y == 0]) + sum(p[y == 1])) / n
pcp_null <- (sum(1 - p0[y0 == 0]) + sum(p0[y0 == 1])) / n0

# percentage correct predictions full model
pcp_full
#> [1] 0.8374718

# percentage correct predictions null model
pcp_null
#> [1] 0.6614221

Создан в 2019-05-29 пакетом Представить (v0.3.0)

Наконец, если вы действительно хотите что-то вроде 0 и 1, вы можете использовать simulate(), который возвращает значения в том же масштабе, что и исходный ответ.Чем вы можете сравнить, сколько из смоделированных ответов соответствует фактическому ответу: 

library(glmmTMB)

data(mtcars)

m <- glmmTMB(formula = vs ~ hp + wt + (1 | gear), family = binomial, data = mtcars)

# simulate response, first column = successes
s <- simulate(m)$sim_1[, 1]

# proportion of response values that equal simulated responses
mean(mtcars$vs == s)
#> [1] 0.875

Создано в 2019-05-29 пакетом представлением (v0.3.0)

0 голосов
/ 23 мая 2019

Функция прогнозирования дает вам вероятность того, что оффшорная переменная является «оффшорной», учитывая предикторы, которые вы указали в модели для значений, выраженных в наборе тестовых данных.Возможно, вы захотите посмотреть на ответы на этот вопрос: Введите параметр функции Forect () .Где они подчеркивают разницу между «Предсказанием.глм» и «Предсказанием.рпарт».Я не могу повторить ваш анализ на основе предоставленных вами данных, но посмотрите на простой пример из: https://www.theanalysisfactor.com/r-tutorial-glm1/. Это показывает, как работает glm (в целом) и что дает вам функция pregnet.glm.Надеюсь, это поможет.

...