Implémenter un test de Vuong

[Véronique Darmendrail]

Bonjour

J'ai changé de PC et version de R ou rstudio et mon ancien script qui utilisait le package "games" avec la fonction Vuong ne peut plus marcher car la package n'existe pas pour cette version de R. Je compare deux régressions linéaires non-nested faites avec lm.
J'ai essayé une alternative. J'ai trouvé le package "nonnest2" avec la fonction testvuong. Mon problème c'est que j'ai a priori des résultats divergents. Que font-ils de différent ou comment puis je mettre les bonnes options pour retrouver la même chose ??? ça m'inquiète un peu.
Je ne vois pas trop comment donner du code reproductible mais voilà mon ancien code et le nouveau :

Code : Tout sélectionner

 library(games) vuongtest(a, b)

et ça donnait ça comme résultat
##
## Vuong test for non-nested models
##
## Model 1 log-likelihood: -20709
## Model 2 log-likelihood: -20696
## Observations: 4382
## Test statistic: -1.2
##
## Neither model is significantly preferred (p = 0.24)

et maintenant :

Code : Tout sélectionner

library(nonnest2)
vuongtest(a, b)

avec ceci comme résultat :
Variance test
H0: Model 1 and Model 2 are indistinguishable
H1: Model 1 and Model 2 are distinguishable
w2 = 0.011, p = 0.0237

Non-nested likelihood ratio test
H0: Model fits are equal for the focal population
H1A: Model 1 fits better than Model 2
z = -1.789, p = 0.963
H1B: Model 2 fits better than Model 1
z = -1.789, p = 0.03678

Sinon, y a-t-il un remplaçant exact de "games" ?

Merci d'avance pour vos éclaircissements

[Logez Maxime]

Bonjour,

tu peux aussi utiliser la fonction vuong.test du package mpath comme juge de paix.

Cordialement,
Maxime

[Véronique Darmendrail]

Bonjour,

tu peux aussi utiliser la fonction vuong.test du package mpath comme juge de paix.

Cordialement,
Maxime

Merci@Maxime, j'avais voulu le faire mais si je ne me trompes pas je n'ai pas des objets du bon type en ayant utilisé lm pour mes régressions ?

[Logez Maxime]

re,

il te suffit de transformer ton lm en glm, sachant que par défaut la fonction glm (family = gaussian) revient à faire un lm.

Cordialement,
Maxime

[Véronique Darmendrail]

Bonsoir
@Maxime : j'ai essayé mais la version que j'ai pu installer de mpath n'a plus la vuong.test : https://www.rdocumentation.org/packages ... s/0.4-2.22

Sinon j'ai pu récupérer mon ancien ordinateur et tester les deux fonctions Vuong des packages "nonnest2" et "games" sur le mêmes données et même ordinateur. Je confirme ma différence de résultat de cette après-midi.

J'ai finalement testé sur ce même ordinateur les deux packages avec les données de test de MASS et j'ai adapté le code d'exemple de nonnest2 pour avoir des régressions linéaires simples et non des lois de poisson. Certes la conclusion est la même avec les deux tests mais la valeur de la statistique est différente (-4.059 versus - 3.6). J'aimerais comprendre pourquoi car ça explique sans doute la différence de conclusions avec mes données et régressions à moi.

Code : Tout sélectionner

require(MASS)


#avec des reg linéaires

house1 <- lm(Freq ~ Infl + Type + Cont, data=housing)
house2 <- lm(Freq ~ Infl + Sat,  data=housing)


vuongtest(house2, house1)
#avec games
library(games)
vuong(house2, house1)

Merci d'avance pour votre aide

[Véronique Darmendrail]

Bonsoir,

J'ai écrit au mainteneur du package "nonest2". Voici sa réponse : "I had a look at this and believe the nonnest2 answer to be correct. This is because the Vuong test is based on casewise log-likelihood contributions. And if you sum up these contributions, you should be able to obtain the model's log-likelihood as reported by logLik() (which is a base function, so I trust it to be the correct answer). Only the nonnest2 results can do this, as shown below.

Code : Tout sélectionner

ll_games <- games:::indivLogLiks(house1)
ll_nn2 <- llcont(house1)

c(sum(ll_games), sum(ll_nn2), logLik(house1))


Comparing the games code to the nonnest2 code, it seems that they differ in the residual SD that is used to compute the likelihood contributions. games appears to directly use the error SD estimated by the model, whereas nonnest2 changes the denominator (multiplying by (N - p)/N).

While I could be wrong, I have found that agreement with logLik() is a pretty good indicator of what is correct.
"
Quelqu'un a-t-il un avis sur ça ?

J'ai quand même voulu comparer avec ce que propose le package "pscl". Ils n'ont pas non plus les mêmes résultats ( pscl propose 3 stat : raw, AIC-corrected et BIC-corrected mais aucun ne correspond à celui de nonnest2) :

Code : Tout sélectionner

require(MASS)
house1 <- glm(Freq ~ Infl + Type + Cont, family=poisson, data=housing)
house2 <- glm(Freq ~ Infl + Sat, family=poisson, data=housing)
library(nonnest2)
vuongtest(house2, house1)

library(pscl)
vuong(house2, house1)

J'ai aussi voulu comparer 3 régressions linéaires simples entre les packages "games", "nonest2" et "pscl". Et je crois que je bloque pour écrire proprement une régression lm avec glm de sorte à ce que les vuong ne soient pas modifiés par cette transformation. Je ne dois pas bien le faire car les résultats ne sont pas les mêmes avec lm ou glm avec nonnest2 et pscl fait un message d'erreur...

Code : Tout sélectionner

#### Installe + charge le package DAAG
#install.packages("DAAG")
library(DAAG)

### Charge les données
data(houseprices)

### Aperçu des premières lignes
head(houseprices)


reg1 <- lm(sale.price ~  area, data=houseprices)
reg2 <- lm(sale.price ~  bedrooms,  data=houseprices)

library(nonnest2)
nonnest2::vuongtest(reg2,reg1)

library(games)
games::vuong(reg2,reg1)

# Si je les écris avec glm


reg3 <-glm(sale.price ~  area, data=houseprices, family=gaussian(link="identity"))
reg4 <- glm(sale.price ~  bedrooms,  data=houseprices, family=gaussian(link="identity"))



library(nonnest2)
nonnest2::vuongtest(reg4,reg3)



library(pscl)
pscl::vuong(reg4,reg3)

# sans mettre link="identity"


reg3 <-glm(sale.price ~  area, data=houseprices, family=gaussian)
reg4 <- glm(sale.price ~  bedrooms,  data=houseprices, family=gaussian)



library(nonnest2)
vuongtest(reg4,reg3)

library(pscl)
pscl::vuong(reg4,reg3)

Merci d'avance pour votre aide sur ces sujets.