Simulieren Sie lineare Regression mit Heteroskedastizität

9

Ich versuche, einen Datensatz zu simulieren, der mit meinen empirischen Daten übereinstimmt, bin mir jedoch nicht sicher, wie ich die Fehler in den Originaldaten abschätzen soll. Die empirischen Daten beinhalten Heteroskedastizität, aber ich bin nicht daran interessiert, sie weg zu transformieren, sondern ein lineares Modell mit einem Fehlerterm zu verwenden, um Simulationen der empirischen Daten zu reproduzieren.

Nehmen wir zum Beispiel an, ich habe einen empirischen Datensatz und ein Modell:

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)

Mit plot(n,y)bekommen wir folgendes. Geben Sie hier die Bildbeschreibung ein

Wenn ich jedoch versuche, die Daten zu simulieren, simulate(mod)wird die Heteroskedastizität entfernt und nicht vom Modell erfasst.

Ich kann ein verallgemeinertes Modell der kleinsten Quadrate verwenden

VMat <- varFixed(~n)
mod2 = gls(y ~ n, weights = VMat)

Das bietet eine bessere Modellanpassung basierend auf AIC, aber ich weiß nicht, wie ich Daten mithilfe der Ausgabe simulieren soll.

Meine Frage ist, wie ich ein Modell erstelle, mit dem ich Daten simulieren kann, die mit den ursprünglichen empirischen Daten übereinstimmen (n und y oben). Insbesondere brauche ich eine Möglichkeit, Sigma2, den Fehler, mithilfe eines Modells abzuschätzen.

user44796
quelle
1
Das lineare Modell erfasst also keine bedingte Heteroskedastizität, es sei denn, es versucht dies explizit mit einem der wenigen Ansätze. Ökonometrische Standardtechniken passen die Standardfehler an Parametern an, um die Heteroskedastizität zu berücksichtigen, modellieren sie jedoch nicht explizit.
generic_user
Du hast recht. Ich versuche, ein lineares Modell zu verwenden, um die Heterogenität zu erfassen. Ich denke, ich sollte ein verallgemeinertes Modell der kleinsten Quadrate verwenden. Wenn es andere Empfehlungen gibt, werde ich sie versuchen.
user44796
Es gibt einen Fehler in Ihrem Code, Sie müssen `lm (y ~ n)`
kjetil b halvorsen
1
Ich verstehe Ihre Frage nicht, weil Ihr Code genau das erreicht, was Sie in seinem Titel zu verlangen scheinen: Er simuliert eine lineare Regression mit heteroskedastischen Fehlern. Fragen Sie nach Methoden, um ein Modell für die Heteroskedastizität abzuschätzen? Wenn ja, müssen Sie ein Modell angeben!
whuber
Hoffentlich habe ich meine Frage mit Änderungen geklärt. In der obigen Frage repräsentieren n und y die empirischen Daten. Ich möchte ein Modell an die Daten anpassen und dann das Modell verwenden, um simulierte Daten zu generieren, die dem Mittelwert und den Residuen der Originaldaten entsprechen.
user44796

Antworten:

9

Um Daten mit einer variierenden Fehlervarianz zu simulieren, müssen Sie den Datengenerierungsprozess für die Fehlervarianz angeben. Wie in den Kommentaren erwähnt, haben Sie dies getan, als Sie Ihre Originaldaten generiert haben. Wenn Sie echte Daten haben und dies versuchen möchten, müssen Sie nur die Funktion identifizieren, die angibt, wie die Restvarianz von Ihren Kovariaten abhängt. Die Standardmethode hierfür besteht darin, Ihr Modell anzupassen, zu überprüfen, ob es angemessen ist (abgesehen von der Heteroskedastizität), und die Residuen zu speichern. Diese Residuen werden zur Y-Variablen eines neuen Modells. Unten habe ich das für Ihren Datengenerierungsprozess getan. (Ich sehe nicht, wo Sie den zufälligen Startwert festgelegt haben, daher sind dies nicht buchstäblich dieselben Daten, sondern sollten ähnlich sein, und Sie können meinen genau reproduzieren, indem Sie meinen Startwert verwenden.)

set.seed(568)  # this makes the example exactly reproducible

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)
res    = residuals(mod)

windows()
  layout(matrix(1:2, nrow=2))
  plot(n,y)
  abline(coef(mod), col="red")
  plot(mod, which=3)

Geben Sie hier die Bildbeschreibung ein

Beachten Sie, dass R‚s ? Plot.lm geben Ihnen einen Plot (vgl hier ) der Quadratwurzel der absoluten Werte der Residuen, helfend mit einem Lowess fit überlagert, was genau das , was Sie brauchen. (Wenn Sie mehrere Kovariaten haben, möchten Sie dies möglicherweise für jede Kovariate separat bewerten.) Es gibt den geringsten Hinweis auf eine Kurve, aber dies sieht so aus, als ob eine gerade Linie die Daten gut anpasst. Passen wir also explizit dieses Modell an:

res.mod = lm(sqrt(abs(res))~fitted(mod))
summary(res.mod)
# Call:
# lm(formula = sqrt(abs(res)) ~ fitted(mod))
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -3.3912 -0.7640  0.0794  0.8764  3.2726 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 1.669571   0.181361   9.206  < 2e-16 ***
# fitted(mod) 0.023558   0.003157   7.461 2.64e-12 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.285 on 198 degrees of freedom
# Multiple R-squared:  0.2195,  Adjusted R-squared:  0.2155 
# F-statistic: 55.67 on 1 and 198 DF,  p-value: 2.641e-12
windows()
  layout(matrix(1:4, nrow=2, ncol=2, byrow=TRUE))
  plot(res.mod, which=1)
  plot(res.mod, which=2)
  plot(res.mod, which=3)
  plot(res.mod, which=5)

Geben Sie hier die Bildbeschreibung ein

Wir brauchen uns keine Sorgen zu machen, dass die Restvarianz auch für dieses Modell im Skalenortungsdiagramm zuzunehmen scheint - das muss im Wesentlichen geschehen. Es gibt wieder den geringsten Hinweis auf eine Kurve, sodass wir versuchen können, einen quadratischen Term anzupassen und zu sehen, ob dies hilft (aber nicht):

res.mod2 = lm(sqrt(abs(res))~poly(fitted(mod), 2))
summary(res.mod2)
# output omitted
anova(res.mod, res.mod2)
# Analysis of Variance Table
# 
# Model 1: sqrt(abs(res)) ~ fitted(mod)
# Model 2: sqrt(abs(res)) ~ poly(fitted(mod), 2)
#   Res.Df    RSS Df Sum of Sq     F Pr(>F)
# 1    198 326.87                          
# 2    197 326.85  1  0.011564 0.007 0.9336

Wenn wir damit zufrieden sind, können wir diesen Prozess jetzt als Add-On verwenden, um Daten zu simulieren.

set.seed(4396)  # this makes the example exactly reproducible
x = n
expected.y = coef(mod)[1] + coef(mod)[2]*x
sim.errors = rnorm(length(x), mean=0,
                   sd=(coef(res.mod)[1] + coef(res.mod)[2]*expected.y)^2)
observed.y = expected.y + sim.errors

Beachten Sie, dass bei diesem Prozess nicht mehr garantiert wird, dass er den tatsächlichen Datengenerierungsprozess findet als bei jeder anderen statistischen Methode. Sie haben eine nichtlineare Funktion verwendet, um die Fehler-SDs zu generieren, und wir haben sie mit einer linearen Funktion approximiert. Wenn Sie den tatsächlichen Datengenerierungsprozess a priori kennen (wie in diesem Fall, weil Sie die Originaldaten simuliert haben), können Sie ihn auch verwenden. Sie können entscheiden, ob die Annäherung hier für Ihre Zwecke gut genug ist. Wir kennen jedoch normalerweise den tatsächlichen Prozess der Datengenerierung nicht und verwenden basierend auf Occams Rasiermesser die einfachste Funktion, die den Daten, die wir angesichts der verfügbaren Informationsmenge angegeben haben, angemessen entspricht. Sie können auch Splines oder schickere Ansätze ausprobieren, wenn Sie dies bevorzugen. Die bivariaten Verteilungen sehen mir ziemlich ähnlich,

Geben Sie hier die Bildbeschreibung ein

gung - Monica wieder einsetzen
quelle
Dies war eigentlich eine Schlussfolgerung, zu der ich anfing zu kommen, aber niemals zu einer so eleganten Antwort gekommen wäre.
user44796
5

Sie müssen die Heteroskedastizität modellieren. Ein Ansatz ist das R-Paket (CRAN) dglm, ein dispersionsverallgemeinertes lineares Modell. Dies ist eine Erweiterung von glm's, die zusätzlich zu den üblichen glmglm ein zweites glm zur Dispersion aus den Resten des ersten glm passt. Ich habe keine Erfahrung mit solchen Modellen, aber sie scheinen vielversprechend ... Hier ist ein Code:

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)

library(dglm)  ### double glm's

mod2   <-  dglm(y ~ n, ~ n, gaussian,ykeep=TRUE,xkeep=TRUE,zkeep=TRUE)
### This uses log link for the dispersion part, should also try identity link ..

y2 <-  simulate(mod2)

plot(n, y2$sim_1)

mod3  <-  dglm(y ~ n, ~ n, gaussian, dlink="identity", ykeep=TRUE,xkeep=TRUE,zkeep=TRUE)  ### This do not work because it leads to negative weights!

Das simulierte Diagramm ist unten dargestellt:

Geben Sie hier die Bildbeschreibung ein

Das Diagramm sieht so aus, als hätte die Simulation die geschätzte Varianz verwendet, aber ich bin mir nicht sicher, da die Funktion simulate () keine Methoden für dglms ...

(Eine andere Möglichkeit, dies zu untersuchen, ist die Verwendung des RPakets gamlss, das einen anderen Ansatz zur Modellierung der Varianz als Funktion von Kovariablen verwendet.)

kjetil b halvorsen
quelle
1
Das doppelt verallgemeinerte lineare Modell scheint die Originaldaten angemessen zu modellieren. Mir ist unklar, wie der Restfehler mit Predict () modelliert wird. Ich werde das untersuchen müssen.
user44796