. use http://www.farys.org/daten/bmi.dta, clear
(ALLBUS 2008)
. reg bmi alter sex einkommen bjahre
Source | SS df MS Number of obs = 1238
-------------+------------------------------ F( 4, 1233) = 34.76
Model | 2719.0338 4 679.75845 Prob > F = 0.0000
Residual | 24114.6464 1233 19.5577018 R-squared = 0.1013
-------------+------------------------------ Adj R-squared = 0.0984
Total | 26833.6802 1237 21.6925466 Root MSE = 4.4224
------------------------------------------------------------------------------
bmi | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
alter | .0536144 .007522 7.13 0.000 .0388571 .0683717
sex | -1.666415 .2686091 -6.20 0.000 -2.193396 -1.139433
einkommen | -.000323 .0001467 -2.20 0.028 -.0006108 -.0000352
bjahre | -.4016991 .0850305 -4.72 0.000 -.5685195 -.2348787
_cons | 30.26551 1.090151 27.76 0.000 28.12675 32.40427
------------------------------------------------------------------------------
. cprplot alter, msize(0.3) jitter(2) lowess lsopts(bwidth(0.4))
. graph export cprplot.png, replace
(file cprplot.png written in PNG format)
Für die Splines scheinen nach graphischer Analyse und einigem Herumprobieren z.B. 41 und 57 Jahre gut als Knotenpunkte zu funktionieren:
. capture drop alter0 alter41 alter57
. gen alter0 = cond(alter<=41, alter, 41) if alter<.
(12 missing values generated)
. gen alter41 = cond(alter<=41, 0, cond(alter<=57, alter, 57) - 41) if alter<.
(12 missing values generated)
. gen alter57 = cond(alter<=57, 0, alter-57) if alter<.
(12 missing values generated)
. reg bmi sex bjahre einkommen alter0 alter41 alter57
Source | SS df MS Number of obs = 1238
-------------+------------------------------ F( 6, 1231) = 34.18
Model | 3831.99638 6 638.666063 Prob > F = 0.0000
Residual | 23001.6838 1231 18.6853646 R-squared = 0.1428
-------------+------------------------------ Adj R-squared = 0.1386
Total | 26833.6802 1237 21.6925466 Root MSE = 4.3227
------------------------------------------------------------------------------
bmi | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
sex | -1.725687 .2633621 -6.55 0.000 -2.242375 -1.208999
bjahre | -.3935012 .0831995 -4.73 0.000 -.5567297 -.2302727
einkommen | -.0004829 .0001482 -3.26 0.001 -.0007736 -.0001922
alter0 | .0777626 .025597 3.04 0.002 .027544 .1279811
alter41 | .1805961 .0278122 6.49 0.000 .1260316 .2351606
alter57 | -.0988631 .0213213 -4.64 0.000 -.1406933 -.057033
_cons | 29.19182 1.29287 22.58 0.000 26.65535 31.72829
------------------------------------------------------------------------------
Man sieht, dass der BMI mit dem Alter zunimmt. Diese Zunahme beschleunigt sich weiter in den Jahren 41-57 bevor sich der Trend deutlich umkehrt. Mit mkspline
können wir die Splines noch etwas bequemer erstellen:
mkspline alter0 41 alter41 57 alter57 = alter
Für die Splines sind natürlich auch andere Knotenpunkte denkbar. V.a. ist die Wahl der Anzahl Knoten arbiträr. Man könnte noch ein wenig herumspielen mit anderen Kombinationen, z.B. indem man sich ein kleines Programm schreibt:
capture program drop myspline
program define myspline
syntax [, knot1(integer 1) knot2(integer 1)]
gen alter0 = cond(alter<=`knot1', alter, `knot1') if alter<.
gen alter`knot1' = cond(alter<=`knot1', 0, cond(alter<=`knot2', alter, `knot2') - `knot1') if alter<.
gen alter`knot2' = cond(alter<=`knot2', 0, alter-`knot2') if alter<.
reg bmi einkommen bjahre sex alter0 alter`knot1' alter`knot2'
capture drop alter0 alter`knot1' alter`knot2'
end
myspline, knot1(41) knot2(57)
Das was wir modelliert haben sieht grafisch etwa so aus:
predict yhat if e(sample)
scatter yhat alter if inrange(alter, 18, 97), msize(0.3) || lfit yhat alter if inrange(alter,18,40) || lfit yhat
> alter if inrange(alter,40,56) || lfit yhat alter if inrange(alter,56,97), name(splines, replace)
graph export splines.png, replace
psline
(ssc install pspline
) und schätzen Sie ein Penalized Spline Modell für das Alter. Kontrollieren Sie wieder Geschlecht, Bildungsjahre und Einkommen. Was ist der Vorteil/Nachteil der automatischen Glättung gegenüber der manuellen Modellierung des Zusammenhangs?
. pspline bmi alter sex einkommen bjahre, noisily
(pilot goodness-of-fit chi2(18) = 72.92; p = 0.0000)
(using penalized model ...)
Performing EM optimization:
Performing gradient-based optimization:
Iteration 0: log restricted-likelihood = -3587.5299
Iteration 1: log restricted-likelihood = -3587.497
Iteration 2: log restricted-likelihood = -3587.497
Computing standard errors:
Mixed-effects REML regression Number of obs = 1238
Group variable: _all Number of groups = 1
Obs per group: min = 1238
avg = 1238.0
max = 1238
Wald chi2(4) = 76.14
Log restricted-likelihood = -3587.497 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
bmi | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
alter | .082938 .1068098 0.78 0.437 -.1264052 .2922813
sex | -1.722625 .2645057 -6.51 0.000 -2.241047 -1.204203
einkommen | -.0004813 .0001485 -3.24 0.001 -.0007722 -.0001903
bjahre | -.3919452 .0833645 -4.70 0.000 -.5553366 -.2285538
_cons | 29.02193 2.51085 11.56 0.000 24.10075 33.9431
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
_all: Identity |
sd(__000004..__00000L)(1) | .092022 .0497018 .0319265 .2652357
-----------------------------+------------------------------------------------
sd(Residual) | 4.324826 .0874726 4.156736 4.499712
------------------------------------------------------------------------------
LR test vs. linear regression: chibar2(01) = 44.63 Prob >= chibar2 = 0.0000
(1) __000004 __000005 __000006 __000007 __000008 __000009 __00000A __00000B __00000C __00000D __00000E __00000F
__00000G __00000H __00000I __00000J __00000K __00000L
. graph export pspline.png, replace
(file pspline.png written in PNG format)
Der manuelle Modellbau hat den Vorteil, dass die einzelnen Teile noch griffig interpretierbar sind (z.b. Steigung innerhalb bestimmter Segmente). Der Vorteil von psplines auf der anderen Seite ist, dass Vorhersagen, v.a. an den Knotenpunkten besser funktionieren. Beispiel alter=57:
. qui reg bmi sex bjahre einkommen alter0 alter41 alter57
. margins, at(alter0==41 alter41==16 alter57=0)
Predictive margins Number of obs = 1238
Model VCE : OLS
Expression : Linear prediction, predict()
at : alter0 = 41
alter41 = 16
alter57 = 0
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons | 28.03416 .2780193 100.84 0.000 27.48925 28.57906
------------------------------------------------------------------------------
. qui pspline bmi alter sex einkommen bjahre
. margins, at(alter=57)
Predictive margins Number of obs = 1238
Expression : Linear prediction, fixed portion, predict()
at : alter = 57
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons | 26.53652 3.771104 7.04 0.000 19.14529 33.92775
------------------------------------------------------------------------------
Das Modell mit linearen Splines überschätzt den BMI v.a. an direkt an den Knotenpunkten. Ein anderer Vorteil kann sein, dass man eine Variable lediglich kontrollieren will ohne sich speziell für deren Einfluss zu interessieren. In dem Fall stellt das Verfahren eine parameter-sparsame Möglichkeit dar, beispielsweise Alter als potentielle Drittvariable zu kontrollieren.