Как получить предсказания модели с помощью функции svyolr - PullRequest
Новая
       18

Как получить предсказания модели с помощью функции svyolr

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

Я использую функцию svyolr из пакета survey для запуска модели порядковой регрессии.Я хотел бы проанализировать модель, проверив прогнозы, которые она сделает для данных моего опроса и для гипотетических наборов данных.

Поскольку, похоже, нет метода predict для svyolrфункция, как я могу получить предсказания модели из подобранной svyolr модели?

Вот некоторые примеры данных и простой пример модели для работы: 

library(survey)

# Load example data
data(api)

# Create a survey design object
dclus1 <- svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)

# Add an ordinal variable to the data
dclus1 <- update(dclus1, mealcat=cut(meals,c(0,25,50,75,100)))

# Fit an ordinal regression model
m <- svyolr(formula = mealcat ~ avg.ed + stype,
            design = dclus1)

# Generate fake data to use in evaluating model predictions
fake_data <- dclus1$variables[,c("avg.ed", "stype")]
fake_data <- subset(fake_data, !is.na(avg.ed))

1 Ответ

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

svyolr основан на MASS :: polr, так что я думаю, вы могли бы взломать алгоритм прогнозирования для этого ..

я сомневаюсь, что доверительные интервалы, которые генерирует glm.predict :: polr.predict (), являются законными 

library(glm.predict)
library(survey)


data = MASS::survey
data$Smoke = ordered(MASS::survey$Smoke,levels=c("Never","Occas","Regul","Heavy"))
model1 = MASS::polr(Smoke ~ Sex + Height, data=data)
summary(model1)
# comparing a man to a woman
polr.predict(model1, c(1,170),sim.count=10000)


this_design <- svydesign( ~ 1 , data = data , weights = ~ Age )
model2 = svyolr(Smoke ~ Sex + Height, this_design)
class(model2) <- 'polr'

# prediction for male with height 170
( this_prediction <- polr.predict(model2, c(1,170),sim.count=10000) )


fake_data <- data[ c( 'Sex' , 'Height' ) ]
fake_data <- fake_data[ complete.cases( fake_data ) , ]
fake_data[ , ] <- sapply( fake_data[ , ] , as.numeric )

# lowering sim.count just so it runs faster
res <- apply( fake_data , 1 , function( w ) polr.predict(model2, w ,sim.count=100) )

# prediction for the first record
res[1:4,1]
# lower CI for the first record
res[5:8,1]
# upper CI for the first record
res[9:12,1]

# prediction of heavy smoking for every record
fake_data[ , 'predicted_heavy_smoking' ] <- res[ 4 , ]
...