J'essaie de simuler un ensemble de données qui correspond aux données empiriques dont je dispose, mais je ne sais pas comment estimer les erreurs dans les données d'origine. Les données empiriques incluent l'hétéroscédasticité, mais je ne suis pas intéressé à les transformer, mais plutôt à utiliser un modèle linéaire avec un terme d'erreur pour reproduire des simulations des données empiriques.
Par exemple, supposons que j'en ai un ensemble de données empiriques et un modèle:
n=rep(1:100,2)
a=0
b = 1
sigma2 = n^1.3
eps = rnorm(n,mean=0,sd=sqrt(sigma2))
y=a+b*n + eps
mod <- lm(y ~ n)
en utilisant plot(n,y)
nous obtenons ce qui suit.
Cependant, si j'essaie de simuler les données, simulate(mod)
l'hétéroscédasticité est supprimée et non capturée par le modèle.
Je peux utiliser un modèle des moindres carrés généralisés
VMat <- varFixed(~n)
mod2 = gls(y ~ n, weights = VMat)
qui fournit un meilleur ajustement du modèle basé sur AIC, mais je ne sais pas comment simuler des données en utilisant la sortie.
Ma question est de savoir comment créer un modèle qui me permettra de simuler des données correspondant aux données empiriques originales (n et y ci-dessus). Plus précisément, j'ai besoin d'un moyen d'estimer sigma2, l'erreur, en utilisant soit un modèle?
la source
Réponses:
Pour simuler des données avec une variance d'erreur variable, vous devez spécifier le processus de génération de données pour la variance d'erreur. Comme cela a été souligné dans les commentaires, vous l'avez fait lorsque vous avez généré vos données d'origine. Si vous avez des données réelles et que vous voulez essayer, il vous suffit d'identifier la fonction qui spécifie comment la variance résiduelle dépend de vos covariables. La façon standard de le faire est d'adapter votre modèle, de vérifier qu'il est raisonnable (autre que l'hétéroscédasticité) et d'économiser les résidus. Ces résidus deviennent la variable Y d'un nouveau modèle. Ci-dessous, je l'ai fait pour votre processus de génération de données. (Je ne vois pas où vous placez la graine aléatoire, donc ce ne seront pas littéralement les mêmes données, mais devraient être similaires, et vous pouvez reproduire la mienne exactement en utilisant ma graine.)
Notez que
R
le ? Plot.lm vous donnera un tracé (cf. ici ) de la racine carrée des valeurs absolues des résidus, superposée utilement avec un ajustement plus bas, ce qui est exactement ce dont vous avez besoin. (Si vous avez plusieurs covariables, vous voudrez peut-être les évaluer séparément pour chaque covariable.) Il y a le moindre indice d'une courbe, mais cela ressemble à une ligne droite qui fait un bon travail d'ajustement des données. Adaptons donc explicitement ce modèle:Nous n'avons pas à nous inquiéter du fait que la variance résiduelle semble également augmenter dans le graphique de localisation de l'échelle pour ce modèle - cela doit essentiellement se produire. Il y a à nouveau le moindre indice d'une courbe, nous pouvons donc essayer d'ajuster un terme au carré et voir si cela aide (mais cela ne fonctionne pas):
Si nous en sommes satisfaits, nous pouvons désormais utiliser ce processus comme module complémentaire pour simuler des données.
Notez que ce processus n'est pas plus garanti pour trouver le véritable processus de génération de données que toute autre méthode statistique. Vous avez utilisé une fonction non linéaire pour générer les SD d'erreur, et nous l'avons approximée avec une fonction linéaire. Si vous connaissez réellement le véritable processus de génération de données a priori (comme dans ce cas, parce que vous avez simulé les données d'origine), vous pourriez aussi bien l'utiliser. Vous pouvez décider si l'approximation ici est suffisante pour vos besoins. Cependant, nous ne connaissons généralement pas le véritable processus de génération de données et, sur la base du rasoir d'Occam, nous utilisons la fonction la plus simple qui correspond adéquatement aux données que nous avons fournies, la quantité d'informations disponibles. Vous pouvez également essayer des splines ou des approches plus sophistiquées si vous préférez. Les distributions bivariées me semblent raisonnablement similaires,
la source
Vous devez modéliser l'hétéroskédasticité. Une approche est via le package R (CRAN)
dglm
, modèle linéaire généralisé de dispersion. Il s'agit d'une extension de glm qui, en plus de l'habituelglm
, s'adapte à un deuxième glm pour la dispersion des résidus du premier glm. Je n'ai aucune expérience avec de tels modèles, mais ils semblent prometteurs ... Voici du code:Le tracé simulé est illustré ci-dessous:
L'intrigue ressemble à la simulation qui a utilisé la variance estimée, mais je ne suis pas sûr, car la fonction simulate () n'a pas de méthodes pour dglm ...
(Une autre possibilité à examiner est d'utiliser le
R
packagegamlss
, qui utilise une autre approche pour modéliser la variance en fonction des covariables.)la source