J'utilise ChoiceModelR pour le logit multinomial hiérarchique. Je veux obtenir des estimations pour l'utilité du bien extérieur (qui suit une distribution normale). Le bien extérieur n'a pas de covariables comme les biens intérieurs - par ex. il ne peut pas avoir de prix ou de marque fictive -, donc j'ai mis la liste (none = TRUE) et n'ajoute pas ce non-choix aux données X (comme décrit dans la documentation de ChoiceModelR) mais uniquement aux données y (choice).ChoiceModelR, none = VRAI, "Erreur dans betadraw [bon,] = newbeta [bon,]"
Les itérations commencent normalement, puis à un moment donné, il arrête et dit
"Error in betadraw[good, ] = newbeta[good, ] : NAs are not allowed in subscripted assignments".
Cela se produit probablement parce que la ligne 388 de la fonction « choicemodelr », la « bonne » est NA indice.
Je regardais quelques questions sur choicemodelr (this, this, this), et aussi sur nas dans l'indice (this, this), mais je pense que mon problème est spécifique à cette fonction dans le sens que probablement quelques entrées dans l'itération, devenez si grand/petit que «bon» deviendra NA.
Voici un exemple très simple. Je génère des données avec 3 produits avec différents attributs attribués. Dans la moitié des périodes, le produit 3 n'est pas offert. Les consommateurs de 2000 ont des préférences - distribuées normalement - sur trois attributs (et une préférence pour le bien extérieur). Erreur Logit ajoutée pour être cohérent avec le modèle. Le bien extérieur est indexé en tant que produit 4 (à la fois lorsque 3 et 2 produits étaient dans l'ensemble de choix).
Comment éviter l'erreur NA? Est-ce que je fais quelque chose de mal, ou est-ce un bug général dans la fonction?
J'ai également cherché des exemples en ligne en mettant l'option none = TRUE, mais je n'ai trouvé aucun reproductible. Peut-être que cette option est seulement la chose problématique car il n'y a pas de problème à récupérer les vrais paramètres si je définis none = FALSE, et je ne laisse pas les clients choisir l'option extérieure.
Ainsi, le code qui se traduit par le bogue NA est la suivante:
library("ChoiceModelR")
library("MASS")
set.seed(36)
# Set demand pars
beta_mu = c(-3,4,1)
beta_sigma = diag(c(1,1,1))
alfa_mu = 5 #outside good mean utility
alfa_sigma = 2 #outside good sd
# Three/two products, 3 vars (2 continuous,1 dummy)
threeprod <- list()
twoprod <- list()
purchase <- list()
for (t in 1:1000){
threeprod[[t]] = cbind(rep(t,3),c(1,1,1),c(1,2,3),runif(3),runif(3),ceiling(runif(3,-0.5,0.5)))
purchase[[t]] = which.max(rbind(threeprod[[t]][,c(4,5,6)]%*%mvrnorm(1,beta_mu,beta_sigma) +
matrix(-log(-log(runif(3))), 3, 1),rnorm(1,alfa_mu,alfa_sigma)))
threeprod[[t]] = cbind(threeprod[[t]],c(purchase[[t]],0,0))
}
for (t in 1001:2000){
twoprod[[t]] = cbind(rep(t,2),c(1,1),c(1,2),runif(2),runif(2),ceiling(runif(2,-0.5,0.5)))
purchase[[t]] = which.max(rbind(twoprod[[t]][,c(4,5,6)]%*%mvrnorm(1,beta_mu,beta_sigma) +
matrix(-log(-log(runif(2))), 2, 1),rnorm(1,alfa_mu,alfa_sigma)))
if (purchase[[t]] == 3) {purchase[[t]] <- 4}
twoprod[[t]] = cbind(twoprod[[t]],c(purchase[[t]],0))
}
X <- rbind(do.call(rbind,threeprod),do.call(rbind,twoprod))
xcoding <- c(1,1,1)
mcmc = list(R = 5000, use = 2000)
options = list(none=TRUE, save=TRUE, keep=5)
out = choicemodelr(X, xcoding, mcmc = mcmc,options = options)