El enlace del registro GAM no funciona sin valores de partida

Question

Estoy tratando de estimar un modelo de regresión GAM utilizando la implementación de gam de la mgcv paquete. Tengo un modelo gaussiano que funciona para la dispersión y un enlace logarítmico para los predictores lineales, pero recibo el error

>"Error in eval(expr, envir, enclos) : cannot find valid starting values: please specify some". 

Edición 1 - La sintaxis exacta es

splineWAR <- gam(WAR ~ s(zAge, bs="cr") + s(zAdjProd, bs="cr") + s(zSOPct, bs="cr") + s(zBBPct, bs="cr"), family=gaussian(link="log"), data = mydata,  start=c(0, 0, 0, 0, 0))

He leído los hilos pertinentes aquí y aquí pero no han podido aplicar los pasos sugeridos a una regresión múltiple. Por ejemplo, cuando intento establecer valores iniciales para las 5 variables de mi regresión (1 dependiente y 4 independientes) añadiendo el start=c(n1, n2, n3, n4, n5) (donde el n son la media de la variable correspondiente), recibo el mismo error aunque aparentemente estoy copiando la sintaxis exactamente del primer enlace. ¿Puede alguien sugerirme qué debo hacer a continuación? Gracias.

Edición 2 El código en el gam.fit que se ejecuta justo antes del error es -

if (!(validmu(mu) && valideta(eta))) 

stop("Can't find valid starting values: please specify some")

Answer 1

Actualización Simon Wood ha corregido este error en mgcv en la versión 1.7-25. La entrada dice:

* bugs fixed whereby etastart etc were not passed to initial.spg and
  get.null.coefs.

---

Hay varios lugares en los que se puede producir un mensaje de error como el que muestras. El primero es en la inicialización de la GAM:

> gaussian()$initialize
    expression({
        n <- rep.int(1, nobs)
        if (is.null(etastart) && is.null(start) && is.null(mustart) && 
            ((family$link == "inverse" && any(y == 0)) || (family$link == 
            "log" && any(y <= 0)))) 
        stop("cannot find valid starting values: please specify some")
    mustart <- y
})

Fíjate en la última cláusula: (family$link == "log" && any(y <= 0))

La primera parte de la cláusula es TRUE en su caso, ¿qué pasa con la segunda parte? Eso fallará, pero la pregunta entonces es, ¿por qué falló como is.null(start) debería haber sido FALSE en su caso. Este código se llama en realidad en gam.fit a través de mgcv:::estimate.gam y por lo tanto mgcv:::initial.spg como muestra el traceback() :

> traceback()
6: stop("cannot find valid starting values: please specify some")
5: eval(expr, envir, enclos)
4: eval(family$initialize)
    3: initial.spg(G$X, G$y, G$w, G$family, G$S, G$off, G$L, G$lsp0)
2: estimate.gam(G, method, optimizer, control, in.out, scale, gamma, 
       ...)
1: gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat, family = gaussian(link = "log"), 
       start = c(1, 2, 3, 4, 5))

Si miramos en mgcv:::initial.spg observamos estas líneas (ignorar el nobs uno)

start <- etastart <- mustart <- NULL
nobs <- nrow(X)
eval(family$initialize)

Es decir, el código anterior evalúa la expresión que he mostrado antes. Pero lo hace después de borrar start .

Creo que se trata de un error ya que no veo cómo se puede encajar el modelo que se quiere con la forma en que está codificado. Por ejemplo, utilizando

 library(mgcv)
 set.seed(2) ## simulate some data... 
 dat <- gamSim(1,n=400,dist="normal",scale=2)
 b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3), data=dat, start = 1:5)

observamos que eso falla, pero nos indica cuántos valores iniciales hay que proporcionar.

> b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3), data=dat, start = 1:5)
Error in gam.fit(G, family = G$family, control = control, gamma = gamma,  : 
  Length of start should equal 37 and correspond to initial coefs.

OK. Cambia esto por un "log" enlazar y reajustar con 37 valores de partida

> b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,
+          family = gaussian(link = "log"), start = runif(37))
Error in eval(expr, envir, enclos) : 
  cannot find valid starting values: please specify some

y sigue fallando.

Una vez que el autor haya solucionado este problema, todavía tendrá que especificar los valores iniciales correctos para el número de términos en su modelo spline.

Answer 2

Un ejemplo de su sintaxis sería útil.

Los GLM requieren tanto una familia (o lo que yo llamo una familia de varianza) como una función de enlace para definir el algoritmo de puntuación de Fisher que resuelve sus estimaciones de parámetros. Con la varianza de Poisson y el enlace logarítmico, se trata de una regresión o regresión de Poisson, que se consigue con el argumento family=poisson a glm . Sin embargo, utilice el siguiente argumento family=binomial(link="log") y se obtiene la regresión del riesgo relativo.

La mayoría de las familias de la R La función GLM permite especificar link="log" como argumento opcional del objeto familia (por ejemplo gaussian , gamma , poisson ). Los GLM irregulares, especialmente aquellos para los que el rango de la función de enlace es mayor que el de la media ajustada en la varianza, tienen extrañas constricciones impuestas en el espacio de los parámetros que la puntuación de Fisher no puede acomodar.

Utilizando el traceback() es siempre útil con errores como estos. También puede encontrar la iteración donde el algoritmo diverge especificando glm.control=list(maxit=1) para un estimador de 1 paso, glm.control=list(maxit=2) para un paso de 2, y así sucesivamente. Trazar su $\beta^{(i)}$ Las estimaciones para la (i)-ésima iteración le ayudarán a ver lo que ocurre antes de que el hessiano se vuelva singular, la puntuación de Fisher diverja y R explote.

Su problema con los valores iniciales puede deberse a que está suministrando las medias de las variables de respuesta cuando las variables son contrasta entre el Registros de los medios. Por lo tanto, los cocientes de los logaritmos de las medias serían un mejor punto de partida. Personalmente, si esto fuera un problema, ajustaría un GLM normal para asegurarme de que el algoritmo no es universalmente divergente y utilizaría las estimaciones de los parámetros de ese modelo para empezar otro.

Por ejemplo: Alimentar las estimaciones de la regresión logística (ORs) para obtener las razones de riesgo (RRs)

## retrospective incidence of something nasty
data <- data.frame(cases=rpois(10, 10), controls=rpois(10, 1000), age=factor(seq(10), labels=c('0-10', '10-20', '20-30', '30-40', '40-50', '50-60', '60-70', '70-80', '80-90', '90-100')))

## logistic regression
fit <- glm(cbind(cases, controls) ~ age, data=data, family=binomial)

## relative risk regression
fit2 <- glm(cbind(cases, controls) ~ age, data=data, family=binomial(link='log'), start=coef(fit))

Answer 3

El parámetro "start" toma valores para los parámetros no para las variables de la regresión. El modelo sólo tiene 4 parámetros (uno para cada variable dependiente) por lo que debe probar start=c(0,0,0,0).