Korelacja między estymatorami OLS dla punktu przecięcia i nachylenia

W prostym modelu regresji

y = β_{0} + β_{1} x + ε,

$y = \beta_0 + \beta_1 x + \varepsilon,$

estymatory OLS i są skorelowane. $\hat{\beta}_0^{OLS}$ $\hat{\beta}_1^{OLS}$

Wzór na korelację między dwoma estymatorami jest (jeśli poprawnie ją wyprowadziłem):

Corr ({\hat{β}}_{0}^{O L S}, {\hat{β}}_{1}^{O L S}) = \frac{- \sum_{i = 1}^{n} x_{i}}{\sqrt{n} \sqrt{\sum_{i = 1}^{n} x_{i}^{2}}} .

$\operatorname{Corr}(\hat{\beta}_0^{OLS},\hat{\beta}_1^{OLS}) = \frac{-\sum_{i=1}^{n}x_i}{\sqrt{n} \sqrt{\sum_{i=1}^{n}x_i^2} }.$

Pytania:

Jakie jest intuicyjne wytłumaczenie obecności korelacji?
Czy obecność korelacji ma jakieś ważne implikacje?

Wpis został zredagowany, a twierdzenie, że korelacja zanika wraz z rozmiarem próbki, zostało usunięte. (Dzięki @whuber i @ChristophHanck.)

regression least-squares estimators

— Richard Hardy
źródło

Wzór jest poprawny, ale czy mógłbyś wyjaśnić, jakich asymptotyków używasz? W końcu w wielu przypadkach korelacja nie zanika - stabilizuje się. Rozważmy np. Eksperyment, w którym są binarne i załóżmy, że dane są gromadzone przez naprzemienne między a . Wtedy a korelacja będzie zawsze zbliżona do , bez względu na to, jak duże będzie .

x_{i}

$x_i$

x_{i}

$x_i$

1

$1$

0

$0$

\sum x_{i} = \sum x_{i}^{2} \approx n / 2

$\sum x_i = \sum x_i^2 \approx n/2$

\sqrt{2} / 2 \neq 0

$\sqrt{2}/2 \ne 0$

n

$n$

— whuber

Powiedziałbym, że znika tylko, jeśli : wpisz co oznacza

E (X) = 0

$E(X)=0$

Corr ({\hat{β}}_{0}^{O L S}, {\hat{β}}_{1}^{O L S}) = \frac{- \frac{1}{N} \sum_{i = 1}^{N} x_{i}}{\sqrt{\frac{N \sum_{i = 1}^{N} x_{i}^{2}}{N^{2}}}} = \frac{- \frac{1}{N} \sum_{i = 1}^{N} x_{i}}{\sqrt{\frac{\sum_{i = 1}^{N} x_{i}^{2}}{N}}},

$\operatorname{Corr}(\hat{\beta}_0^{OLS},\hat{\beta}_1^{OLS}) = \frac{-\frac{1}{N}\sum_{i=1}^{N}x_i}{\sqrt{\frac{N\sum_{i=1}^{N}x_i^2}{N^2}}} = \frac{-\frac{1}{N}\sum_{i=1}^{N}x_i}{\sqrt{\frac{\sum_{i=1}^{N}x_i^2}{N}}},$

- E (X) / \sqrt{E (X^{2})}

$-E(X)/\sqrt{E(X^2)}$

— Christoph Hanck

Rzeczywiście, pominąłem jeden gdy wyprowadzałem zachowanie korelacji wraz ze wzrostem . Więc Whuber i ChristophHanck mają rację. Nadal interesuje mnie intuicyjne wyjaśnienie, dlaczego korelacja jest niezerowa, i wszelkie użyteczne implikacje . (Ja nie powiedzieć, że korelacja powinna intuicyjnie wynosić zero, po prostu nie mają intuicję tutaj.)

n

$n$

n

$n$

— Richard Hardy

Twoja formuła starannie pokazuje np., Że w przypadku regresora o wartości średniej korelacja z punktem przecięcia zanika.

x

$x$

— Michael M

Powiązane: Dlaczego standardowy błąd przechwytywania zwiększa się o dalsze od zera?

\bar{x}

$\bar x$

— Gung - Przywróć Monikę

Pozwól mi spróbować w następujący sposób (naprawdę nie jestem pewien, czy to przydatna intuicja):

Na podstawie mojego powyższego komentarza korelacja będzie z grubsza wynosić

- \frac{E (X)}{\sqrt{E (X^{2})}}

$-\frac{E(X)}{\sqrt{E(X^2)}}$ Zatem, jeśli

E (X) > 0

$E(X)>0$ zamiast

E (X) = 0

$E(X)=0$ , większość danych zostanie skupiona w prawo od zera. Zatem jeśli współczynnik nachylenia staje się większy, formuła korelacji stwierdza, że punkt przecięcia musi być mniejszy - co ma pewien sens.

Myślę o czymś takim:

W niebieskiej próbce oszacowanie nachylenia jest bardziej płaskie, co oznacza, że ocena przecięcia może być większa. Nachylenie złotej próbki jest nieco większe, więc punkt przecięcia może być nieco mniejszy, aby to zrekompensować.

Z drugiej strony, jeśli , możemy mieć dowolne nachylenie bez żadnych ograniczeń na przecięciu. $E(X)=0$

Mianownik formuły można również interpretować według tych zasad: jeśli dla danej średniej zmienność mierzona za pomocą wzrasta, dane są rozmazane na osi , tak że skutecznie „wygląda” więcej średnią zero znowu poluzowanie ograniczeń na przecięcia na dany średnią . $E(X^2)$ $x$ $X$

Oto kod, który, mam nadzieję, w pełni wyjaśnia rysunek:

n <- 30
x_1 <- sort(runif(n,2,3))
beta <- 2
y_1 <- x_1*beta + rnorm(n) # the golden sample

x_2 <- sort(runif(n,2,3)) 
beta <- 2
y_2 <- x_2*beta + rnorm(n) # the blue sample

xax <- seq(-1,3,by=.001)
plot(x_1,y_1,xlim=c(-1,3),ylim=c(-4,7),pch=19,col="gold",ylab="y",xlab="x")
abline(lm(y_1~x_1),col="gold",lwd=2)
abline(v=0,lty=2)
lines(xax,beta*xax) # the "true" regression line
abline(lm(y_2~x_2),col="lightblue",lwd=2)
points(x_2,y_2,pch=19,col="lightblue")

— Christoph Hanck
źródło

Aby uzyskać praktyczną implikację, rozważ opracowanie i zastosowanie krzywej kalibracyjnej dla instrumentu laboratoryjnego. W celu opracowania kalibracji znane wartości

są testowane za pomocą przyrządu i mierzone są wartości wyjściowe przyrządu

, a następnie regresja liniowa. Następnie na przyrząd nakładana jest nieznana próbka, a nowa wartość

jest używana do przewidywania nieznanego

na podstawie kalibracji regresji liniowej. Analiza błędu oszacowania nieznanego

wymagałaby korelacji między oszacowaniami nachylenia regresji i przecięcia.

x

$x$

y

$y$

y

$y$

x

$x$

x

$x$

— EdM

Możesz podążać za wprowadzeniem Dougherty'ego do ekonometrii , być może biorąc pod uwagę, że jest zmienną niestochastyczną, i określając średnie odchylenie kwadratowe jako $x$ $x$ . Należy zauważyć, że MSD jest mierzony w kwadracie jednostek(np. Jeślijest wto MSD jest w), podczas gdy średnie odchylenie kwadratowe pierwiastka, $\DeclareMathOperator{\MSD}{MSD}\MSD(x) = \frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})^2$ $x$ $x$ $\text{cm}$ $\text{cm}^2$ jest w oryginalnej skali. To daje $\DeclareMathOperator{\RMSD}{RMSD}\RMSD(x)=\sqrt{\MSD(x)}$

Corr ({\hat{β}}_{0}^{O L S}, {\hat{β}}_{1}^{O L S}) = \frac{- \bar{x}}{\sqrt{MSD (x) + {\bar{x}}^{2}}}

$\DeclareMathOperator{\Corr}{Corr}\Corr(\hat{\beta}_0^{OLS},\hat{\beta}_1^{OLS}) = \frac{-\bar{x}}{\sqrt{\MSD(x) + \bar{x}^2}}$

To powinno pomóc zobaczyć jak korelacja ma wpływ zarówno na średnią z (w szczególności, korelacja między nachylenie i przecięcie estymatorów jest usuwany, jeśli zmienna jest wyśrodkowany), a także przez jego rozprzestrzeniania . (Ten rozkład mógł również uczynić asymptotyki bardziej oczywistymi!) $x$ $x$

Ponownie podkreślę znaczenie tego wyniku: jeśli nie ma oznaczającego zera, możemy go przekształcić, odejmując , aby był teraz wyśrodkowany. Jeśli dopasujemy linię regresji na szacunki nachylenia i przecięcia są nieskorelowane - niedoszacowanie lub przeszacowanie w jednym nie powoduje zwykle niedoszacowania lub przeszacowania w drugim. Ale ta linia regresji jest po prostu tłumaczeniem linii regresji na ! Błąd standardowy punktu przecięcia z na linii jest po prostu miarą niepewności $x$ $\bar{x}$ $y$ $x - \bar{x}$ $y$ $x$ $y$ $x - \bar{x}$ $\hat y$ kiedy twoja przetłumaczona zmienna ; gdy ta linia jest tłumaczony z powrotem do jego pierwotnego położenia, ten powraca do stanu błąd standardowy przy . Bardziej ogólnie, błąd standardowy w każdej wartości jest tylko błąd standardowy przecięcia regresji w odpowiednio tłumaczone ; Błąd standardowy przy jest oczywiście błędu standardowego punktu przecięcia w oryginalnych, nie podlegającym translacji regresji. $x - \bar x = 0$ $\hat y$ $x = \bar x$ $\hat y$ $x$ $y$ $x$ $\hat y$ $x=0$

Ponieważ możemy tłumaczyć , w pewnym sensie nie ma nic szczególnego w , a zatem nic szczególnego . Przy odrobinie myślenia, co mam do powiedzenia prace dla na każdej wartości , co jest przydatne, jeśli szukasz wgląd np przedziały ufności dla średnich odpowiedzi od linii regresji. Jednak widzieliśmy, że tam jest coś szczególnego przy , bo jest tu, że błędy w szacowanej wysokości linii regresji - co jest oczywiście szacuje się na $x$ $x=0$ $\hat \beta_0$ $\hat y$ $x$ $\hat y$ $x=\bar x$ - a błędy w szacowanym nachyleniu linii regresji nie mają ze sobą nic wspólnego. Szacowany przechwytujący jest i błędy w estymacji musi wynikać albo z oszacowania lub estymacji(ponieważ uważanyjak wolny stochastyczny); teraz wiemy, że te dwa źródła błędów są nieskorelowane, algebraicznie jasne jest, dlaczego powinna istnieć ujemna korelacja między szacowanym nachyleniem a przecięciem (przeszacowanie nachylenia będzie miało tendencję do niedoceniania przechwytywania, o ile $\bar y$ $\hat \beta_0 = \bar y - \hat \beta_1 \bar x$ $\bar y$ $\hat \beta_1$ $x$ ), ale dodatnia korelacja między szacowaną przecięcia i estymowanej odpowiedzi przy . Ale widzą takie relacje również bez algebry. $\bar x < 0$ $\hat y = \bar y$ $x = \bar x$

Wyobraź sobie oszacowaną linię regresji jako linijkę. Ten władca musi przejść . Właśnie widzieliśmy, że istnieją dwie zasadniczo niezwiązane ze sobą niepewności w lokalizacji tej linii, które wizualizuję kinestetycznie jako niepewność „drgającą” i „równoległą przesuwającą się”. Zanim otoczysz linijkę, przytrzymaj ją w $(\bar x, \bar y)$ $(\bar x, \bar y)$ jako czop, a następnie nadaj mu mocny odgłos związany z twoją niepewnością na zboczu. Linijka będzie się dobrze kołysać, bardziej gwałtownie, więc jeśli jesteś bardzo niepewny co do nachylenia (rzeczywiście, poprzednio dodatnie nachylenie będzie całkiem prawdopodobne, że będzie ujemne, jeśli twoja niepewność jest duża), ale zwróć uwagę, że wysokość linii regresji przy pozostaje niezmieniony przez tego rodzaju niepewność, a efekt twangu jest bardziej zauważalny w dalszej odległości od środka, na który patrzysz. $x=\bar x$

Aby „zsunąć” linijkę, chwyć ją mocno i przesuń w górę i w dół, uważając, aby równolegle do pierwotnej pozycji - nie zmieniać nachylenia! To, jak energicznie przesuwać go w górę i w dół, zależy od tego, jak niepewna jest wysokość linii regresji, gdy przechodzi ona przez punkt średni; pomyśl o tym, jaki byłby standardowy błąd punktu przecięcia, gdyby został przetłumaczony, aby oś przechodziła przez punkt średni. Alternatywnie, ponieważ szacowana wysokość linii regresji wynosi po prostu , jest to również błąd standardowy . Zauważ, że ten rodzaj „przesuwającej się” niepewności wpływa na wszystkie punkty linii regresji w jednakowy sposób, w przeciwieństwie do „twangu”. $x$ $y$ $\bar y$ $\bar y$

Te dwie niewiadome zastosowanie niezależnie (dobrze, uncorrelatedly, ale jeśli założymy, zazwyczaj dystrybuowane warunki błędach to powinny one być technicznie niezależny) tak Wysokości wszystkich punktów na linii regresji są dotknięte przez „twanging” niepewności, który jest zero u wredne i coraz gorzej, a „przesuwająca się” niepewność, która jest wszędzie taka sama. (Czy widzisz związek z przedziałami ufności regresji, które obiecałem wcześniej, szczególnie to, że ich szerokość jest węższa przy ?) $\hat y$ $\bar x$

Obejmuje to niepewność przy , co jest w zasadzie to, co rozumiemy przez standardowego błędu . Załóżmy teraz, że jest po prawej stronie ; następnie zmiana wykresu na wyższy szacowany spadek ma tendencję do zmniejszania naszego szacowanego punktu przecięcia, co ujawni szybki szkic. Jest to korelacja ujemna przewidywana przez $\hat y$ $x=0$ $\hat \beta_0$ $\bar x$ $x=0$ gdyjest dodatnie. I odwrotnie, jeślijest na lewo od, zobaczysz, że wyższe oszacowane nachylenie ma tendencję do zwiększania naszego szacowanego przecięcia, zgodnie zdodatniąkorelacją, którą twoje równanie przewiduje, gdyjest ujemne. Zauważ, że jeślijest daleko od zera, ekstrapolacja linii regresji niepewnego gradientu w kierunku $\frac{-\bar{x}}{\sqrt{\MSD(x) + \bar{x}^2}}$ $\bar x$ $\bar x$ $x=0$ $\bar x$ $\bar x$ $y$ -osi stają się coraz bardziej niepewne (amplituda „twang” pogarsza się od średniej). Przycisków „twanging” błąd w perspektywie znacznie przewyższają „ślizgową” błąd w perspektywie więc błąd jest prawie całkowicie określona przez każdy błąd . Jak można łatwo zweryfikować algebraicznie, jeśli weźmiemy bez zmiany MSD lub odchylenie standardowe błędy , korelacja pomiędzy $- \hat \beta_1 \bar x$ $\bar y$ $\hat \beta_0$ $\hat \beta_1$ $\bar x \to \pm \infty$ $s_u$ $\hat \beta_0$ ma tendencję do. $\hat \beta_1$ $\mp 1$

Aby to zilustrować (możesz kliknąć obraz prawym przyciskiem myszy i zapisać go lub wyświetlić go w pełnym rozmiarze w nowej karcie, jeśli ta opcja jest dla Ciebie dostępna). Wybrałem rozważenie powtarzania próbek , gdzie to iid, w stosunku do ustalonego zestawu wartości przy , więc $y_i = 5 + 2x_i + u_i$ $u_i \sim N(0, 10^2)$ $x$ $\bar x = 10$ $\mathbb{E}(\bar y)=25$ . W tym układzie istnieje dość silna korelacja ujemna między szacowanym nachyleniem i przecięciem, a słabsza dodatnia korelacja między , szacunkową średnią odpowiedzią przy i oszacowanym przecięciem. Animacja pokazuje kilka symulowanych próbek, z linią regresji próbki (złota) narysowaną nad rzeczywistą (czarną) linią regresji. Drugi rząd pokazuje, jak wyglądałby zbiór oszacowanych linii regresji, gdyby wystąpił błąd tylko w szacowanym a nachylenia pasowały do prawdziwego nachylenia (błąd „przesuwu”); wtedy, jeśli wystąpił błąd tylko na zboczach i $\bar y$ $x=\bar x$ $\bar y$ $\bar y$ dopasował jego wartość populacji (błąd „wykręcania”); i wreszcie, jak rzeczywiście wyglądał zbiór linii szacunkowych, gdy oba źródła błędów zostały połączone. Zostały one zakodowane kolorami według wielkości faktycznie oszacowanego przechwytywania (nie przechwytywania pokazanego na pierwszych dwóch wykresach, na których wyeliminowano jedno ze źródeł błędów) od niebieskiego dla niskich przechwyceń do czerwonego dla wysokich przechwyceń. Należy pamiętać, że od samych kolorach widzimy, że próbki z niskiego tendencję do wytwarzania obniżyć szacowane przechwytuje, podobnie jak próbki z wysokiego $\bar y$ szacowane stoki. Następny rząd pokazuje symulowane (histogram) i teoretyczne (krzywa normalna) rozkłady prób oszacowań, a ostatni rząd pokazuje wykresy rozrzutu między nimi. Zauważ, że nie ma korelacji między a szacowanym nachyleniem, ujemna korelacja między szacowanym przecięciem a nachyleniem oraz dodatnia korelacja między przecięciem a . $\bar y$ $\bar y$

Co robi MSD w mianowniku ? Rozkładanie zakreswartości zmierzyć się jest dobrze znany, co pozwala na bardziej precyzyjne oszacowanie nachylenia, a intuicja jest jasne od szkicu, ale nie pozwalają oszacowaćlepiej. Sugeruję wizualizację zbliżenia MSD do zera (tj. Punkty próbkowania tylko bardzo blisko średniej), aby twoja niepewność na zboczu stała się ogromna: pomyśl wielkie wielkie twangi, ale bez zmiany niepewności poślizgu. Jeśli twojaośjest jakąkolwiek odległością od(innymi słowy, jeśli $\frac{-\bar{x}}{\sqrt{\MSD(x) + \bar{x}^2}}$ $x$ $\bar y$ $x$ $y$ $\bar x$ $\bar x \neq 0$ ) przekonasz się, że niepewność w twoim punkcie przechwytywania całkowicie zdominowana jest przez błąd skrętu związany ze spadkiem. W przeciwieństwie do tego, jeśli zwiększysz rozpiętość swoich pomiarów , bez zmiany średniej, znacznie poprawisz precyzję oszacowania nachylenia i potrzebujesz tylko delikatnych twangów do linii. Wysokość przechwytywania jest teraz zdominowana przez niepewność poślizgu, która nie ma nic wspólnego z szacowanym nachyleniem. Jest to zgodne z faktem algebraicznym, że korelacja między szacowanym nachyleniem i przecięciem dąży do zera jako a gdy , w kierunku $x$ $\MSD(x) \to \pm \infty$ $\bar x \neq 0$ $\pm 1$ (znak jest przeciwieństwem znaku ) jako . $\bar x$ $\MSD(x) \to 0$

Korelacja estymatorów nachylenia i przecięcia była funkcją zarówno i MSD (lub RMSD) , więc w jaki sposób zwiększają się ich względne udziały? W rzeczywistości liczy się tylko stosunek do RMSD . Zgodnie z intuicją geometryczną RMSD daje nam rodzaj „jednostki naturalnej” dla ; jeśli przeskalujemy oś za pomocą to jest to odcinek poziomy, który pozostawia szacowany punkt przecięcia i bez zmian, daje nam nowy $\bar x$ $x$ $\bar x$ $x$ $x$ $x$ $w_i = x_i / \RMSD(x)$ $\bar y$ i mnoży szacunkowe nachylenie przez RMSD . Wzór na korelację między nowymi estymatorami nachylenia i przecięcia dotyczy tylko , który wynosi jeden, i , który jest stosunkiem $\RMSD(w)=1$ $x$ $\RMSD(w)$ $\bar w$ . Ponieważ oszacowanie przecięcia pozostało niezmienione, a oszacowanie nachylenia zostało jedynie pomnożone przez stałą dodatnią, korelacja między nimi nie uległa zmianie: stąd korelacja międzypierwotnymnachyleniem a przecięciem musi również zależeć tylko od $\frac{\bar x}{\RMSD(x)}$ $\frac{\bar x}{\RMSD(x)}$ $\frac{-\bar x}{\sqrt{\MSD(x)+\bar{x}^2}}$ $\RMSD(x)$ $\Corr\left(\hat \beta_0, \hat \beta_1 \right) = \frac{- (\bar x / \RMSD(x))}{\sqrt{1 + (\bar x / \RMSD(x))^2}}$ .

To find the correlation between $\hat \beta_0$ and $\bar y$ , consider $\DeclareMathOperator{\Cov}{Cov}\Cov(\hat \beta_0, \bar y)=\Cov(\bar y - \hat \beta_1 \bar x, \bar y)$ . By bilinearity of $\Cov$ this is $\Cov(\bar y, \bar y) - \bar x \Cov(\hat \beta_1, \bar y)$ . The first term is $\operatorname{Var}(\bar y)=\frac{\sigma_u^2}{n}$ while the second term we established earlier to be zero. From this we deduce

Corr ({\hat{β}}_{0}, \bar{y}) = \frac{1}{\sqrt{1 + (\bar{x} / RMSD (x))^{2}}}

$\Corr(\hat \beta_0, \bar y)=\frac{1}{\sqrt{1 + (\bar x/\RMSD(x))^2}}$

So this correlation also depends only on the ratio $\frac{\bar x}{\RMSD(x)}$ . Note that the squares of $\Corr(\hat \beta_0, \hat \beta_1)$ and $\Corr(\hat \beta_0, \bar y)$ sum to one: we expect this since all sampling variation (for fixed $x$ ) in $\hat \beta_0$ is due either to variation in $\hat \beta_1$ or to variation in $\bar y$ , and these sources of variation are uncorrelated with each other. Here is a plot of the correlations against the ratio $\frac{\bar x}{\RMSD(x)}$ .

The plot clearly shows how when $\bar x$ is high relative to the RMSD, errors in the intercept estimate are largely due to errors in the slope estimate and the two are closely correlated, whereas when $\bar x$ is low relative to the RMSD, it is error in the estimation of $\bar y$ that predominates, and the relationship between intercept and slope is weaker. Note that the correlation of intercept with slope is an odd function of the ratio $\frac{\bar x}{\RMSD(x)}$ , so its sign depends on the sign of $\bar x$ and it is zero if $\bar x=0$ , whereas the correlation of intercept with $\bar y$ is always positive and is an even function of the ratio, i.e. it doesn't matter what side of the $y$ -axis that $\bar x$ is. The correlations are equal in magnitude if $\bar x$ is one RMSD away from the $y$ -axis, when $\Corr(\hat \beta_0, \bar y)=\frac{1}{\sqrt{2}} \approx 0.707$ and $\Corr(\hat \beta_0, \hat \beta_1)=\pm \frac{1}{\sqrt{2}} \approx \pm 0.707$ where the sign is opposite that of $\bar x$ . In the example in the simulation above, $\bar x=10$ and $\RMSD(x) \approx 5.16$ so the mean was about $1.93$ RMSDs from the $y$ -axis; at this ratio, the correlation between intercept and slope is stronger, but the correlation between intercept and $\bar y$ is still not negligible.

As an aside, I like to think of the formula for the standard error of the intercept,

s . e . ({\hat{β}}_{0}^{O L S}) = \sqrt{s_{u}^{2} (\frac{1}{n} + \frac{{\bar{x}}^{2}}{n MSD (x)})}

$\operatorname{s.e.}(\hat \beta_0^{OLS}) = \sqrt{s_u^2 \left( \frac{1}{n} + \frac{{\bar x}^2 }{n \MSD(x)} \right) }$

as $\sqrt{\text{sliding error} + \text{twanging error}}$ , and ditto for the formula for the standard error of $\hat y$ at $x = x_0$ (used for confidence intervals for the mean response, and of which the intercept is just a special case as I explained earlier via a translation argument),

s . e . (\hat{y}) = \sqrt{s_{u}^{2} (\frac{1}{n} + \frac{(x_{0} - \bar{x})^{2}}{n MSD (x)})}

$\operatorname{s.e.}(\hat y) = \sqrt{s_u^2 \left( \frac{1}{n} + \frac{(x_0 - \bar x)^2}{n \MSD(x)} \right) }$

R code for plots

require(graphics)
require(grDevices)
require(animation

#This saves a GIF so you may want to change your working directory
#setwd("~/YOURDIRECTORY")
#animation package requires ImageMagick or GraphicsMagick on computer
#See: http://www.inside-r.org/packages/cran/animation/docs/im.convert
#You might only want to run up to the "STATIC PLOTS" section
#The static plot does not save a file, so need to change directory.

#Change as desired
simulations <- 100 #how many samples to draw and regress on
xvalues <- c(2,4,6,8,10,12,14,16,18) #used in all regressions
su <- 10 #standard deviation of error term
beta0 <- 5 #true intercept
beta1 <- 2 #true slope
plotAlpha <- 1/5 #transparency setting for charts
interceptPalette <- colorRampPalette(c(rgb(0,0,1,plotAlpha),
            rgb(1,0,0,plotAlpha)), alpha = TRUE)(100) #intercept color range
animationFrames <- 20 #how many samples to include in animation

#Consequences of previous choices
n <- length(xvalues) #sample size
meanX <- mean(xvalues) #same for all regressions
msdX <- sum((xvalues - meanX)^2)/n #Mean Square Deviation
minX <- min(xvalues)
maxX <- max(xvalues)
animationFrames <- min(simulations, animationFrames)

#Theoretical properties of estimators
expectedMeanY <- beta0 + beta1 * meanX
sdMeanY <- su / sqrt(n) #standard deviation of mean of Y (i.e. Y hat at mean x)
sdSlope <- sqrt(su^2 / (n * msdX))
sdIntercept <- sqrt(su^2 * (1/n + meanX^2 / (n * msdX)))


data.df <- data.frame(regression = rep(1:simulations, each=n),
                      x = rep(xvalues, times = simulations))

data.df$y <- beta0 + beta1*data.df$x + rnorm(n*simulations, mean = 0, sd = su) 

regressionOutput <- function(i){ #i is the index of the regression simulation
  i.df <- data.df[data.df$regression == i,]
  i.lm <- lm(y ~ x, i.df)
  return(c(i, mean(i.df$y), coef(summary(i.lm))["x", "Estimate"],
          coef(summary(i.lm))["(Intercept)", "Estimate"]))
}

estimates.df <- as.data.frame(t(sapply(1:simulations, regressionOutput)))
colnames(estimates.df) <- c("Regression", "MeanY", "Slope", "Intercept")

perc.rank <- function(x) ceiling(100*rank(x)/length(x))
rank.text <- function(x) ifelse(x < 50, paste("bottom", paste0(x, "%")), 
                                paste("top", paste0(101 - x, "%")))
estimates.df$percMeanY <- perc.rank(estimates.df$MeanY)
estimates.df$percSlope <- perc.rank(estimates.df$Slope)
estimates.df$percIntercept <- perc.rank(estimates.df$Intercept)
estimates.df$percTextMeanY <- paste("Mean Y", 
                                    rank.text(estimates.df$percMeanY))
estimates.df$percTextSlope <- paste("Slope",
                                    rank.text(estimates.df$percSlope))
estimates.df$percTextIntercept <- paste("Intercept",
                                    rank.text(estimates.df$percIntercept))

#data frame of extreme points to size plot axes correctly
extremes.df <- data.frame(x = c(min(minX,0), max(maxX,0)),
              y = c(min(beta0, min(data.df$y)), max(beta0, max(data.df$y))))

#STATIC PLOTS ONLY

par(mfrow=c(3,3))

#first draw empty plot to reasonable plot size
with(extremes.df, plot(x,y, type="n", main = "Estimated Mean Y"))
invisible(mapply(function(a,b,c) { abline(a, b, col=c) }, 
                 estimates.df$Intercept, beta1, 
                 interceptPalette[estimates.df$percIntercept]))

with(extremes.df, plot(x,y, type="n", main = "Estimated Slope"))
invisible(mapply(function(a,b,c) { abline(a, b, col=c) }, 
                 expectedMeanY - estimates.df$Slope * meanX, estimates.df$Slope, 
                 interceptPalette[estimates.df$percIntercept]))

with(extremes.df, plot(x,y, type="n", main = "Estimated Intercept"))
invisible(mapply(function(a,b,c) { abline(a, b, col=c) }, 
                 estimates.df$Intercept, estimates.df$Slope, 
                 interceptPalette[estimates.df$percIntercept]))

with(estimates.df, hist(MeanY, freq=FALSE, main = "Histogram of Mean Y",
                        ylim=c(0, 1.3*dnorm(0, mean=0, sd=sdMeanY))))
curve(dnorm(x, mean=expectedMeanY, sd=sdMeanY), lwd=2, add=TRUE)

with(estimates.df, hist(Slope, freq=FALSE, 
                        ylim=c(0, 1.3*dnorm(0, mean=0, sd=sdSlope))))
curve(dnorm(x, mean=beta1, sd=sdSlope), lwd=2, add=TRUE)

with(estimates.df, hist(Intercept, freq=FALSE, 
                        ylim=c(0, 1.3*dnorm(0, mean=0, sd=sdIntercept))))
curve(dnorm(x, mean=beta0, sd=sdIntercept), lwd=2, add=TRUE)

with(estimates.df, plot(MeanY, Slope, pch = 16,  col = rgb(0,0,0,plotAlpha), 
                        main = "Scatter of Slope vs Mean Y"))

with(estimates.df, plot(Slope, Intercept, pch = 16, col = rgb(0,0,0,plotAlpha),
                        main = "Scatter of Intercept vs Slope"))

with(estimates.df, plot(Intercept, MeanY, pch = 16, col = rgb(0,0,0,plotAlpha),
                        main = "Scatter of Mean Y vs Intercept"))


#ANIMATED PLOTS

makeplot <- function(){for (i in 1:animationFrames) {

  par(mfrow=c(4,3))

  iMeanY <- estimates.df$MeanY[i]
  iSlope <- estimates.df$Slope[i]
  iIntercept <- estimates.df$Intercept[i]

  with(extremes.df, plot(x,y, type="n", main = paste("Simulated dataset", i)))
  with(data.df[data.df$regression==i,], points(x,y))
  abline(beta0, beta1, lwd = 2)
  abline(iIntercept, iSlope, lwd = 2, col="gold")

  plot.new()
  title(main = "Parameter Estimates")
  text(x=0.5, y=c(0.9, 0.5, 0.1), labels = c(
    paste("Mean Y =", round(iMeanY, digits = 2), "True =", expectedMeanY),
    paste("Slope =", round(iSlope, digits = 2), "True =", beta1),
    paste("Intercept =", round(iIntercept, digits = 2), "True =", beta0)))

  plot.new()
  title(main = "Percentile Ranks")
  with(estimates.df, text(x=0.5, y=c(0.9, 0.5, 0.1),
                          labels = c(percTextMeanY[i], percTextSlope[i],
                                     percTextIntercept[i])))


  #first draw empty plot to reasonable plot size
  with(extremes.df, plot(x,y, type="n", main = "Estimated Mean Y"))
  invisible(mapply(function(a,b,c) { abline(a, b, col=c) }, 
                   estimates.df$Intercept, beta1, 
                   interceptPalette[estimates.df$percIntercept]))
  abline(iIntercept, beta1, lwd = 2, col="gold")

  with(extremes.df, plot(x,y, type="n", main = "Estimated Slope"))
  invisible(mapply(function(a,b,c) { abline(a, b, col=c) }, 
                expectedMeanY - estimates.df$Slope * meanX, estimates.df$Slope, 
                interceptPalette[estimates.df$percIntercept]))
  abline(expectedMeanY - iSlope * meanX, iSlope,
         lwd = 2, col="gold")

  with(extremes.df, plot(x,y, type="n", main = "Estimated Intercept"))
  invisible(mapply(function(a,b,c) { abline(a, b, col=c) }, 
                   estimates.df$Intercept, estimates.df$Slope, 
                   interceptPalette[estimates.df$percIntercept]))
  abline(iIntercept, iSlope, lwd = 2, col="gold")

  with(estimates.df, hist(MeanY, freq=FALSE, main = "Histogram of Mean Y",
                          ylim=c(0, 1.3*dnorm(0, mean=0, sd=sdMeanY))))
  curve(dnorm(x, mean=expectedMeanY, sd=sdMeanY), lwd=2, add=TRUE)
  lines(x=c(iMeanY, iMeanY),
        y=c(0, dnorm(iMeanY, mean=expectedMeanY, sd=sdMeanY)),
        lwd = 2, col = "gold")

  with(estimates.df, hist(Slope, freq=FALSE, 
                          ylim=c(0, 1.3*dnorm(0, mean=0, sd=sdSlope))))
  curve(dnorm(x, mean=beta1, sd=sdSlope), lwd=2, add=TRUE)
  lines(x=c(iSlope, iSlope), y=c(0, dnorm(iSlope, mean=beta1, sd=sdSlope)),
        lwd = 2, col = "gold")

  with(estimates.df, hist(Intercept, freq=FALSE, 
                          ylim=c(0, 1.3*dnorm(0, mean=0, sd=sdIntercept))))
  curve(dnorm(x, mean=beta0, sd=sdIntercept), lwd=2, add=TRUE)
  lines(x=c(iIntercept, iIntercept),
        y=c(0, dnorm(iIntercept, mean=beta0, sd=sdIntercept)),
        lwd = 2, col = "gold")

  with(estimates.df, plot(MeanY, Slope, pch = 16,  col = rgb(0,0,0,plotAlpha), 
                          main = "Scatter of Slope vs Mean Y"))
  points(x = iMeanY, y = iSlope, pch = 16, col = "gold")

  with(estimates.df, plot(Slope, Intercept, pch = 16, col = rgb(0,0,0,plotAlpha),
                          main = "Scatter of Intercept vs Slope"))
  points(x = iSlope, y = iIntercept, pch = 16, col = "gold")

  with(estimates.df, plot(Intercept, MeanY, pch = 16, col = rgb(0,0,0,plotAlpha),
                          main = "Scatter of Mean Y vs Intercept"))
  points(x = iIntercept, y = iMeanY, pch = 16, col = "gold")

}}

saveGIF(makeplot(), interval = 4, ani.width = 500, ani.height = 600)

For the plot of correlation versus ratio of $\bar x$ to RMSD:

require(ggplot2)

numberOfPoints <- 200
data.df  <- data.frame(
  ratio = rep(seq(from=-10, to=10, length=numberOfPoints), times=2),
  between = rep(c("Slope", "MeanY"), each=numberOfPoints))
data.df$correlation <- with(data.df, ifelse(between=="Slope",
  -ratio/sqrt(1+ratio^2),
  1/sqrt(1+ratio^2)))

ggplot(data.df, aes(x=ratio, y=correlation, group=factor(between),
                    colour=factor(between))) +
  theme_bw() + 
  geom_line(size=1.5) +
  scale_colour_brewer(name="Correlation between", palette="Set1",
                      labels=list(expression(hat(beta[0])*" and "*bar(y)),
                              expression(hat(beta[0])*" and "*hat(beta[1])))) +
  theme(legend.key = element_blank()) +
  ggtitle(expression("Correlation of intercept estimates with slope and "*bar(y))) +
  xlab(expression("Ratio of "*bar(X)/"RMSD(X)")) +
  ylab(expression(paste("Correlation")))

— Silverfish
źródło

The "twang" and "slide" are my terms. This is my own visual intuition, and not one I have ever seen in any textbook, though the basic ideas here are all standard material. Goodness knows if there is a more technical name than "twang" and "slide"! I based this answer, from memory, on an answer to a related question that I never quite got round to finishing and posting. That had more instructive graphs, which (if I can track down the R code on my old computer, or find the time to reproduce) I will add.

— Silverfish

What a job! Thank you very much! Now my understanding must be in much better shape.

— Richard Hardy

@RichardHardy I have put a simulation animation in, which ought to make things a bit clearer.

— Silverfish