Robię symulację regresji liniowej z R.

Rozważam model regresji jest y_i = a + b_1 * x_1i + b_2 * x_2i + e_i.

Projekt parametrów w następujący sposób:

X_1I ~ N (2,1), X_2I ~ Poisson (4), E_I ~ N (0, 1), THETA = (A, B_1, B_2)

Poniższy kod, który robię, jest to, że chciałbym wygenerować 100 niezależnych próbek losowych (Y, X_1, X_2) 1000 razy, używając dystrybucji, o której wspomniałem powyżej, a także chcę oszacować THETA_HAT (Estymator THETA). Po uzyskaniu THETA_HAT chciałbym wykreślić dystrybucję estymatora (A_HAT), odpowiednio b_1 (B_1_HAT), B_2 (B_2_HAT).

## Construct 1000 x_1
x_1_1000 <- as.data.frame(replicate(n = 1000,expr = rnorm(n = 100, 
  mean = 2, sd = 1)))
colnames(x_1_1000) <- paste("x_1", 1:1000, sep = "_")

x_2_1000 <- as.data.frame(replicate(n = 1000,expr = rpois(n = 100, 
  lambda = 4)))
colnames(x_2_1000) <- paste("x_2", 1:1000, sep = "_")

error_1000 <- as.data.frame(replicate(n = 1000, expr = rnorm(n = 100,
  mean = 0, sd = 1)))
colnames(error_1000) <- paste("e", 1:1000, sep = "_")

y_1000 <- as.data.frame(matrix(data = 0, nrow = 100, ncol = 1000))
y_1000 = 1 + x_1_1000 * 1 + x_2_1000*(-2) + error_1000
colnames(y_1000) <- paste("y", 1:1000, sep = "_")
######################################################################
lms <- lapply(1:1000, function(x) lm(y_1000[,x] ~ x_1_1000[,x] + x_2_1000[,x]))
theta_hat_1000 <- as.data.frame(sapply(lms, coef))

Po wykonaniu regresji liniowej, wystarczy przechowywać wynik do LM, który jest listą. Ponieważ chcę tylko dane współczynnika, przechowuję te współczynniki symulacyjne na "theta_hat_1000", gdy chcę wykreślić wykres dystrybucyjny, nie mogę dostać tego, czego chcę w finale. Próbowałem dwa sposoby rozwiązywania problemu, ale nadal mylić.

Pierwszy sposób, w jaki próbowałem, jest to, że zmieniam nazwę ramki danych "THETA_HAT_1000". Z powodzeniem przemianę przemianę na Column_i, gdzie I od 1 do 1000. Jednak po prostu nie mogę pomyślnie zmienić nazwy wierszy.

rownames(theta_hat_1000[1,]) <- "ahat"
rownames(theta_hat_1000[2,]) <- "x1hat"
rownames(theta_hat_1000[3,]) <- "x2hat"

Powyższy kod nie pokazał żadnego komunikatu o błędzie, ale ostatecznie nie udało się zmienić nazw rzędów. Tak więc próbowałem następującego kodu

rownames(theta_hat_1000) <- c("ahat", "x1hat", "x2hat")

To z powodzeniem przemianowano. Jednak kiedy chcę sprawdzić, że w ramce danych raportuje "NULL"

theta_hat_1000$ahat

NULL

Dlatego zauważyłem, że jest coś dziwnego. Tak więc próbowałem drugiego sposobu jak następujące.

Próbowałem umieścić "THETA_HAT_1000", który jest lista przechowywana w moim globalnym środowisku. Jednak po wykonaniu takich rzeczy, nie dostaję tego, czego chcę. Oczekiwany wynik jest po prostu uzyskanie trzech wierszy i każdego wiersza z 1000 wartości, ale rzeczywiste jest to, że dostałem 3000 Obs z 1 kolumną.

Idealnym wynikiem jest uzyskanie trzech kolumn i każdej kolumny z 1000 wartości i wprowadzenie ich do ramy danych, aby wykonać dalszy proces, taki jak Korzystanie z GGLOT do wykazania dystrybucji szacowanych współczynników.

Utknąłem na tym przez kilka godzin. Byłoby docenione, jeśli ktoś może mi pomóc i dać mi pewne sugestie.

0
Ricky.L 29 październik 2020, 16:44

1 odpowiedź

Najlepsza odpowiedź

Ta linia theta_hat_1000$ahat w kodzie nie działa, ponieważ "Co" to nazwa wiersza, a nie nazwa kolumny w ramce danych. Otrzymasz wynik, dzwoniąc theta_hat_1000["ahat",].

Rozumiem jednak, że pożądany wynik jest właściwie dataframe z 3 kolumnami (i 1000 wierszami) reprezentujących 3 parametry modelu regresji (przechwytuje, X1, X2). Ta linia w kodzie as.data.frame(sapply(lms, coef)) wytwarza dataframe za pomocą 3 rzędów i 1000 kolumn. Możesz na przykład, transponować matrycę przed zmianą go w ramce danych, aby uzyskać 1000 wierszy i 3 kolumn.

theta_hat_1000 <- sapply(lms, coef)
theta_hat_1000 <- as.data.frame(t(theta_hat_1000))
colnames(theta_hat_1000) <- c("ahat", "x1hat", "x2hat")

head(theta_hat_1000)
       ahat     x1hat     x2hat
1 2.0259326 0.7417404 -2.111874
2 0.7827929 0.9437324 -1.944320
3 1.1034906 1.0091594 -2.035405
4 0.9677150 0.8168757 -1.905367
5 1.0518646 0.9616123 -1.985357
6 0.8600449 1.0781489 -2.017061

Teraz możesz również zadzwonić do zmiennych z theta_hat_1000$ahat.

1
AndreasM 29 październik 2020, 14:30