Je suis un étudiant travaillant sur un modèle épidémiologique en R, en utilisant des méthodes de maximum de vraisemblance. J'ai créé ma fonction de vraisemblance logarithmique négative. Il est une sorte de brut à la recherche, mais la voici:Comment obtenir des intervalles de confiance sans inverser une matrice de Hesse singulière dans R?
NLLdiff = function(v1, CV1, v2, CV2, st1 = (czI01 - czV01), st2 = (czI02 - czV02), st01 = czI01, st02 = czI02, tt1 = czT01, tt2 = czT02) {
prob1 = (1 + v1 * CV1 * tt1)^(-1/CV1)
prob2 = (1 + v2 * CV2 * tt2)^(-1/CV2)
-(sum(dbinom(st1, st01, prob1, log = T)) + sum(dbinom(st2, st02, prob2, log = T)))
}
La raison pour laquelle la première ligne a l'air si terrible parce que la plupart des données qu'il faut y est entrée. czI01, par exemple, est déjà déclaré. Je l'ai fait simplement pour que mes derniers appels à la fonction n'aient pas tous d'horribles vecteurs en eux. J'ai ensuite optimisé pour CV1, CV2, v1 et v2 en utilisant mle2 (bibliothèque bbmle). C'est aussi un peu brut à la recherche, et ressemble à:
ml.cz.diff = mle2 (NLLdiff, start=list(v1 = vguess, CV1 = cguess, v2 = vguess, CV2 = cguess), method="L-BFGS-B", lower = 0.0001)
Maintenant, tout fonctionne bien jusqu'à ici. ml.cz.diff me donne des valeurs que je peux transformer en intrigue qui correspond raisonnablement à mes données. J'ai aussi plusieurs modèles différents, et je peux obtenir des valeurs AICc pour les comparer. Cependant, lorsque j'essaie d'obtenir des intervalles de confiance autour de v1, CV1, v2 et CV2, j'ai des problèmes. Fondamentalement, j'obtiens une borne négative sur CV1, ce qui est impossible car il représente en fait un nombre carré dans le modèle biologique ainsi que quelques avertissements.
Existe-t-il un meilleur moyen d'obtenir des intervalles de confiance? Ou, vraiment, un moyen d'obtenir des intervalles de confiance qui ont un sens ici? Ce que je vois arriver, c'est que, par coïncidence, ma matrice hessienne est singulière pour certaines valeurs dans l'espace d'optimisation. Mais, puisque j'optimise sur 4 variables et que je n'ai pas une connaissance de programmation trop étendue, je ne peux pas trouver une bonne méthode d'optimisation qui ne repose pas sur le Hessian. J'ai googlé le problème - il a suggéré que mon modèle est mauvais, mais je suis en train de reconstruire un peu de travail avant que cela suggère que mon modèle n'est vraiment pas terrible (les complots que je fais en utilisant le ml.cz.diff ressemblent aux complots du travail original). J'ai également lu les parties pertinentes du manuel ainsi que le livre de Bolker Ecological Models in R. J'ai également essayé différentes méthodes d'optimisation, ce qui a entraîné une durée d'exécution plus longue mais les mêmes erreurs. La méthode "SANN" n'a pas fini de fonctionner dans une heure, donc je n'ai pas attendu pour voir le résultat.
En un mot: mes intervalles de confiance sont mauvais. Existe-t-il un moyen relativement simple de les corriger dans R?
Mes vecteurs sont les suivants:
czT01 = c(5, 5, 5, 5, 5, 5, 5, 25, 25, 25, 25, 25, 25, 25, 50, 50, 50, 50, 50, 50, 50)
czT02 = c(5, 5, 5, 5, 5, 10, 10, 10, 10, 10, 25, 25, 25, 25, 25, 50, 50, 50, 50, 50, 75, 75, 75, 75, 75)
czI01 = c(25, 24, 22, 22, 26, 23, 25, 25, 25, 23, 25, 18, 21, 24, 22, 23, 25, 23, 25, 25, 25)
czI02 = c(13, 16, 5, 18, 16, 13, 17, 22, 13, 15, 15, 22, 12, 12, 13, 13, 11, 19, 21, 13, 21, 18, 16, 15, 11)
czV01 = c(1, 4, 5, 5, 2, 3, 4, 11, 8, 1, 11, 12, 10, 16, 5, 15, 18, 12, 23, 13, 22)
czV02 = c(0, 3, 1, 5, 1, 6, 3, 4, 7, 12, 2, 8, 8, 5, 3, 6, 4, 6, 11, 5, 11, 1, 13, 9, 7)
et je reçois mes suppositions par:
v = -log((c(czI01, czI02) - c(czV01, czV02))/c(czI01, czI02))/c(czT01, czT02)
vguess = mean(v)
cguess = var(v)/vguess^2
Il est également possible que je fais quelque chose d'autre complètement faux, mais mes résultats semblent raisonnables, donc je havre de paix Je l'ai attrapé.