Сравнение моделей GLM с использованием прогноза

Question

Предположим, у меня есть две модели, созданные путем вызова glm() по тем же данным, но с различными формулами и/или семействами. Теперь я хочу сравнить, какая модель лучше, предсказывая по неизвестным данным. Что-то вроде этого:

mod1 <- glm(formula1, family1, data) 
mod2 <- glm(formula2, family2, data) 
mu1 <- predict(mod1, newdata, type = "response") 
mu2 <- predict(mod2, newdata, type = "response") 

1. Как я могу сказать, какой из предсказаний mu1 или mu2 лучше?
2. Есть ли какая-то простая команда для вычисления логарифмической вероятности предсказания?

Answer 1

Было бы проще ответить на этот вопрос reproducible example.

Это часто имеет смысл выбирать семейство a priori, а не соответствовать слишком хорошей подгонке - например, если у вас есть счетные (неотрицательные целые) ответы без очевидной верхней границы, ваш единственный реальный выбор которая лежит строго внутри экспоненциального семейства, - это Пуассон.

set.seed(101) 
x <- runif(1000) 
mu <- exp(1+2*x) 
y <- rgamma(1000,shape=3,scale=mu/3) 
d <- data.frame(x,y) 

Новые данные:

nd <- data.frame(x=runif(100)) 
nd$y <- rgamma(100,shape=3,scale=exp(1+2*nd$x)/3) 

Fit Gamma и гауссовских:

mod1 <- glm(y~x,family=Gamma(link="log"),data=d) 
mod2 <- glm(y~x,family=gaussian(link="log"),data=d) 

Предсказания:

mu1 <- predict(mod1, newdata=nd, type="response") 
mu2 <- predict(mod2, newdata=nd, type="response") 

Извлечение параметров формы/масштаба:

sigma <- sqrt(summary(mod2)$dispersion) 
shape <- MASS::gamma.shape(mod1)$alpha 

Root среднеквадратическая ошибка:

rmse <- function(x1,x2) sqrt(mean((x1-x2)^2)) 
rmse(mu1,nd$y) ## 5.845 
rmse(mu2,nd$y) ## 5.842 

Негативные вероятностями журнала:

-sum(dgamma(nd$y,shape=shape,scale=mu1/shape,log=TRUE)) ## 276.84 
-sum(dnorm(nd$y,mean=mu2,sd=sigma,log=TRUE)) ## 318.4