Mejorar la estimación de los ceros y cómo tratar la concurrencia en el GAM NB

Question

Estoy tratando de modelar el número de muertes para cada edad, para 4 regiones, a través de 5 años. Ajusté varios modelos, pero terminé con el siguiente modelo:

m1 <- gam(deaths ~ ti(region, bs='re') + ti(ageCenter, bs='cr', k=40) + 
        ti(ageCenter, region, k=c(40,4)) + ageCenter*yearCenter + offset(log(PopMedia)), 
        family = nb, data = male, method = "REML")

El modelo es razonable, todas las variables son significativas. Sin embargo, la concurrencia no es grande para la región:

         para ti(Regiao) ti(ageCenter) ti(ageCenter,Regiao)
worst       1  0.8335333   1.000000000         0.0635483181
observed    1  0.8335333   0.007150858         0.0003614262
estimate    1  0.8335333   0.029291067         0.0005498903

Sin embargo, comparé el número de muertes observado con el previsto y la estimación era baja para una región (la más pequeña). Intenté mejorar el modelo y obtuve una estimación más cercana para la región más pequeña, pero otra región disminuyó el número estimado de ceros. Este es el nuevo modelo:

m2 <- gam(Obitosdx ~ ti(Regiao, bs='re') + ti(ageCenter, bs='cr', k=40, by=Regiao) + 
        ti(ageCenter, Regiao, k=c(40,4)) + ageCenter*yearCenter + offset(log(PopMedia)), 
        family = nb, data = male, method = "REML")

Todas las variables siguen siendo significativas, la concurrencia se mantiene sin cambios ya que todavía hay ti(region) en el modelo. Los residuos entre los dos modelos son similares.

Mi primera pregunta se refiere a los modelos. No estoy seguro de si es apropiado tener ti(ageCenter, bs='cr', k=40, by=Regiao) en el mismo modelo que ti(ageCenter, Regiao, k=c(40,4)) . ¿No sería el ti la interacción de los productos también da diferentes alisados para cada nivel de la Región?

Mi segunda pregunta es la concurrencia para ti(region) ¿preocupante? Si la variable es significativa, ¿cómo tratarla?

Estos son los dos modelos, comparando el m1 y m2 para la región más pequeña y una región intermedia. ¿Tiene alguna sugerencia de mejora? La estimación del número de ceros es bastante baja. Para la región más pequeña, el m2 estima un número de ceros más cercano, mientras que para la región intermedia, el m1 estima un número de ceros más cercano. Cualquier sugerencia será bienvenida.

EDITAR

Siguiendo la sugerencia de @Gavin, he recodificado Region a una variable de factor, que debería haber sido un factor para empezar. Sin embargo, sigo obteniendo una gran desviación para las edades jóvenes. ¿Alguna sugerencia? El gráfico de arriba con la región intermedia corresponde a Alentejo, mientras que la más pequeña a RAM.

Answer 1

No se puede ajustar una suavidad a una variable factorial y esto significa que su ti(ageCenter, Regiao, k=c(40,4)) , ti(Regiao, bs = "re") y ti(ageCentre, by = Regiao) no hacen lo que se espera de ellos.

El hecho de que la primera funcione indica que Region (o Regiao ¿Cuál es? ¿Importa?) no se almacena como variable de factor. Eso significa que el by suave es en realidad un efecto lineal de Region modificado por un suave de ageCentre y no suavizados por separado para las regiones.

La alta concurrencia de los términos se debe probablemente a los problemas anteriores.

También vas a tener que ser muy cuidadoso con el ageCentre*yearCentre ya que esto incluirá un efecto lineal de ageCentre y en algún lugar entre las bases de las splines tendrás un efecto lineal de ageCentre también, en algún lugar del espacio nulo de una de las splines. Esto va a introducir algunos problemas de identificabilidad; o puede, a menos que mgcv está siendo muy inteligente para eliminar esos efectos (y no creo que lo esté haciendo).

Suponiendo que Region es un factor (si no lo es, recodifícalo en uno), entonces un modelo con suavizados separados de ageCentre para cada región más el efecto de ageCentre y yearCentre podría escribirse como

m1 <- gam(deaths ~ ti(yearCenter) +
                   region +
                   ti(ageCenter, bs='cr', k=40, by = region) + 
                   ti(ageCenter, yearCenter) +
                   offset(log(PopMedia)), 
          family = nb, data = male, method = "REML")

Hice el ageCentre por yearCentre una cosa suave, y dudo que quieras estimar un efecto aleatorio con cuatro niveles (pero si debes hacerlo, pon el re base y suave de nuevo en el region término).