Por lo tanto, estoy tratando de comparar dos modelos, fit1 y fit2.

Inicialmente, sólo estaba haciendo anova(fit1,fit2), y esto produjo una salida que he entendido (incluyendo un p-valor).

Sin embargo, cuando me cambié mis modelos de lm () basado en modelos glm () basado en modelos, análisis de la varianza(fit1,fit2) ahora se rindieron Residual Grados de Libertad, Residuos Deviances, y en el Df, Deviances, que estoy teniendo problemas en la interpretación (recursos explicar estos indicadores parecen escasos). Tenía la esperanza de extraer un valor de p para la comparación entre los dos modelos, pero por alguna razón anova(fit1,fit2, prueba=’Chisq’) no está funcionando. Alguna sugerencia?

Me doy cuenta de que, dependiendo de la función de enlace en mi glms, Chi-cuadrado puede no ser el más apropiado para la prueba, pero he utilizado ‘F’ en los contextos apropiados, así como con similar decepción.

Es este problema familiar a alguien más? Sugerencias? Muchas gracias!

Ejemplo:

make_and_compare_models <- function(fitness_trait_name, data_frame_name, vector_for_multiple_regression, predictor_for_single_regression, fam){
fit1<-glm(formula=as.formula(paste(fitness_trait_name,"~", paste(vector_for_multiple_regression, sep="+"))), family=fam, data=data_frame_name)
print ("summary fit 1")
print(summary(fit1))
fit2<- glm(data=data_frame_name, formula=as.formula(paste(fitness_trait_name,"~",predictor_for_single_regression)), family=fam)
print("summary fit 2")
print(summary(fit2))
print("model comparison stats:")
mod_test<-anova(fit2,fit1)
##suggestion #1
print(anova(fit2,fit1, test="Chisq"))
#suggestion #2
print ("significance:")
print (1-pchisq( abs(mod_test$Deviance[2]),df=abs(mod_test$Df[2])))
}
data<-structure(list(ID = c(1L, 2L, 4L, 7L, 9L, 10L, 12L, 13L, 14L, 
15L, 16L, 17L, 18L, 20L, 21L, 22L, 23L, 24L, 25L, 27L, 28L, 29L, 
31L, 34L, 37L, 38L, 39L, 40L, 41L, 43L, 44L, 45L, 46L, 47L, 48L, 
49L, 52L, 55L, 56L, 59L, 60L, 61L, 62L, 63L, 65L, 66L, 67L, 68L, 
69L, 71L), QnWeight_initial = c(158L, 165L, 137L, 150L, 153L, 
137L, 158L, 163L, 159L, 151L, 145L, 144L, 157L, 144L, 133L, 148L, 
151L, 151L, 147L, 158L, 178L, 164L, 134L, 151L, 148L, 142L, 127L, 
179L, 162L, 150L, 151L, 153L, 163L, 155L, 163L, 170L, 149L, 165L, 
128L, 134L, 145L, 147L, 148L, 160L, 131L, 155L, 169L, 143L, 123L, 
151L), Survived_eclosion = c(0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 
1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), Days_wrkr_eclosion_minus20 = c(NA, 
1L, NA, 3L, 0L, 2L, 0L, 1L, 0L, 0L, 0L, 1L, NA, 0L, 7L, 1L, 0L, 
1L, 0L, 1L, 2L, 2L, NA, 2L, 3L, 2L, 2L, NA, 0L, 1L, NA, NA, 0L, 
0L, 0L, 0L, 3L, 3L, 3L, 1L, 0L, 2L, NA, 1L, 0L, 1L, 1L, 3L, 1L, 
2L), MLH = c(0.5, 0.666666667, 0.555555556, 0.25, 1, 0.5, 0.333333333, 
0.7, 0.5, 0.7, 0.5, 0.666666667, 0.375, 0.4, 0.5, 0.333333333, 
0.4, 0.375, 0.3, 0.5, 0.3, 0.2, 0.4, 0.875, 0.6, 0.4, 0.222222222, 
0.222222222, 0.6, 0.6, 0.3, 0.4, 0.714285714, 0.4, 0.3, 0.6, 
0.4, 0.7, 0.625, 0.555555556, 0.25, 0.5, 0.5, 0.6, 0.25, 0.428571429, 
0.3, 0.25, 0.375, 0.555555556), Acon5 = c(0.35387674, 0.35387674, 
0.35387674, 0.35387674, 0.35387674, 0.35387674, 0.35387674, 0, 
0, 1, 0, 1, 0.35387674, 0, 0, 0.35387674, 1, 1, 0, 0, 0, 1, 0, 
0.35387674, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 
0, 0, 1, 0, 0, 0, 1, 0, 0.35387674), Baez = c(1, 1, 1, 0.467836257, 
1, 1, 0, 0, 1, 1, 0, 0.467836257, 1, 0, 0, 0, 0, 1, 0, 0, 0, 
0, 0, 0.467836257, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 
1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1), C294 = c(0, 1, 0, 0, 1, 
0.582542694, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 
0, 1, 1, 0, 0, 0.582542694, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1), C316 = c(1, 1, 0, 0, 0.519685039, 
0.519685039, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0.519685039, 0, 
1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0.519685039, 1, 0, 1, 
1, 0, 0.519685039, 1, 0.519685039, 1, 1, 1, 0.519685039, 0.519685039, 
0, 0.519685039, 0.519685039, 0), i_120_PigTail = c(1, 1, 0, 1, 
0.631236443, 0.631236443, 1, 1, 1, 1, 1, 0, 0.631236443, 1, 1, 
1, 0, 0.631236443, 1, 1, 1, 0, 0, 1, 1, 1, 0.631236443, 0, 1, 
1, 0, 1, 0.631236443, 1, 0, 1, 0, 0, 1, 0.631236443, 0.631236443, 
0, 1, 0, 0.631236443, 0.631236443, 1, 0.631236443, 0.631236443, 
1), i129 = c(0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L), Jackstraw_PigTail = c(0L, 1L, 1L, 0L, 
1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 
0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), Neil_Young = c(0.529636711, 
0, 1, 0, 0.529636711, 0.529636711, 1, 1, 0, 1, 1, 1, 0, 0, 1, 
1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 
1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1), Ramble = c(0, 0, 0, 
0, 0.215163934, 0.215163934, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 
0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.215163934, 0, 
0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0.215163934, 0, 0, 0, 0), Sol_18 = c(1, 
0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 
0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0.404669261, 
1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1)), .Names = c("ID", "QnWeight_initial", 
"Survived_eclosion", "Days_wrkr_eclosion_minus20", "MLH", "Acon5", 
"Baez", "C294", "C316", "i_120_PigTail", "i129", "Jackstraw_PigTail", 
"Neil_Young", "Ramble", "Sol_18"), class = "data.frame", row.names = c(NA, 
-50L))
make_and_compare_models("QnWeight_initial", data, c("Acon5","Baez","C294","C316","i_120_PigTail","i129","Jackstraw_PigTail","Neil_Young","Ramble","Sol_18"), "MLH", "gaussian")
anova(fit1,fit2,test="Chisq") debe de trabajo, a menos que los modelos anidados pasar a tener la misma se ajusta. Puedes dar más detalles?
PS no es la función de enlace, pero la familia que determina si usted debe usar Chi-cuadrado o F (específicamente, si el parámetro de escala es fijo [Poisson, binomial] o estimado [Gaussiano, Gamma, cuasi-verosimilitud se adapta]
gracias por la aclaración. Sólo para estar seguro, es de la Chi-cuadrado para la escala fija los parámetros y F para los cálculos? Además, el resultado de anova(fit1,fit2, test=»Chisq») se obtiene un Pr(<Chi) que no es limitada por (0,1). En otras palabras, no tengo idea de cómo interpretar los valores como -18.215 (también hay grandes números positivos). Ojalá pudiera recordar si este era el problema original que tuve con el uso de test=»Chisq», pero ya no puede.
También, hay un test=»F» analógica? No puedo encontrar nada acerca de la prueba como un parámetro para el análisis de la varianza() en el manual de…
Su ejemplo muestra que usted está comparando no-anidada modelos: el df diferencia (que se muestra en la Df columna) es cero!!! Todos los anova() marco (como se explica en las respuestas abajo) se enmarca alrededor de anidados modelos. Si desea comparar la bondad de ajuste de no-modelos anidados, puede utilizar AIC (con precaución) o el test de Vuong …

OriginalEl autor Atticus29 | 2012-11-05

2 Comentarios

  1. 8

    La diferencia en la desviación entre un «más grande» o de modelos más complejos y de un entramado o «reducido» modelo distribuido (asintóticamente) como una chi-cuadrado de la variable aleatoria, con la diferencia en los grados de libertad de los dos modelos. Por lo que se podría extraer de la desviación y la estimación de la diferencia en los grados de libertad y comparar que a pchisq( la desviación, diff(df) ). El «p-value» es sólo 1 menos ese valor.

    > 1-pchisq(3.84,1)
    [1] 0.05004352

    Si ejecuta el primer ejemplo en el glm de la página de ayuda y, a continuación, agregue un modelo reducido sin el «tratamiento» de la variable, se obtiene:

    glm.D93.o <- glm(counts ~ outcome, family=poisson())
    anova.res <-anova(glm.D93, glm.D93.o)
    anova.res
    #------------
    Analysis of Deviance Table
    Model 1: counts ~ outcome + treatment
    Model 2: counts ~ outcome
    Resid. Df Resid. Dev Df    Deviance
    1         4     5.1291               
    2         6     5.1291 -2 -2.6645e-15
    #---------------
    str(anova.res)
    Classes ‘anova’ and 'data.frame':   2 obs. of  4 variables:
    $ Resid. Df : num  4 6
    $ Resid. Dev: num  5.13 5.13
    $ Df        : num  NA -2
    $ Deviance  : num  NA -2.66e-15
    - attr(*, "heading")= chr  "Analysis of Deviance Table\n" "Model 1: counts ~ outcome + treatment\nModel 2: counts ~ outcome"

    Así que después de mirar cómo las cosas fueron almacenados en el objeto en sí, este da el valor de p para «resultado»:

     1-pchisq( abs(anova.res$Deviance[2]), abs(anova.res$Df[2]))
    [1] 1

    Y este sería el correspondiente procedimiento en el tratamiento+resultado del modelo frente al tratamiento-sólo modelo:

    > glm.D93.t <- glm(counts ~ treatment, family=poisson())
    > anova.res2 <-anova(glm.D93, glm.D93.t)
    > 1-pchisq( abs(anova.res2$Deviance[2]), abs(anova.res2$Df[2]))
    [1] 0.06547071
    Gracias, DWin! Eso responde a mi pregunta!
    el 1-pchisq() no puede ser correcta. He de ejecutar simulaciones con completamente revuelto de datos (es decir, no debe haber ninguna diferencia significativa entre los dos modelos, porque ni el modelo predice con éxito la respuesta), y el informe de los p-valor es siempre «0». Está usted seguro de que no sólo pchisq() en este caso?
    Estoy bastante seguro de que 1-pchisq(3.84,1) devuelve 0.05. Usted necesita asegurarse de que usted está poniendo el valor absoluto de la correcta desviación diferencia en el primer argumento y el correcto grados de libertad en el segundo. El orden del modelo argumentos voltear el signo de la anova $Deviance resultados, pero la abs() debe tener cuidado de que.
    Uno de los puntos. Valores absolutos están ahí. Hmm… ok acabo específicamente designado como el segundo argumento como «df=modelo$DF[2]», y que lo aclaró. Interesante…
    Acepto la corrección… todavía no funciona.

    OriginalEl autor 42-

  2. 1

    Si su 2 modelos anidada, entonces usted puede utilizar el cambio en la desviación de los 2 modelos para ver si el modelo que contiene los parámetros adicionales que se obtiene un mejor ajuste. Si el modelo 1 contiene k parámetros y el modelo 2, contiene los mismos k parámetros más un adicional de m parámetros, entonces el cambio en la desviación de la siguiente manera (aproximadamente) chi-cuadrado de distribución con m grados de libertad. Usted puede utilizar este estadístico de prueba para ver si el modelo 2 es una mejora en el modelo 1.

    Si usted es nuevo en esta zona, recomiendo la lectura de un texto introductorio en GLMs

    perfecto, excepto que no estoy muy seguro de cómo aplicar en la práctica. I. e., sabes por casualidad el R sintaxis?
    Por desgracia, hacía años que he utilizado R. que yo recuerde, el glm.resumen de la salida se utiliza para proporcionar todo lo que se necesitaba para este cálculo. Esperamos que usted conseguirá un R específicos de la respuesta, en vez de sólo teórica.

    OriginalEl autor mathematician1975

Dejar respuesta

Please enter your comment!
Please enter your name here