Statistische modellen: overzicht
Help! Statistiek! Doorlopende serie laagdrempelige lezingen, voor iedereen vrij toegankelijk.
- eenvoudig voorbeeld (nut van een model)
Doel: Informeren over statistiek in klinisch onderzoek. Tijd:
- enkele eenvoudige regressie-modellen - lineaire regressie - logit - probit - kwantiel-regressie in R
Derde woensdag in de maand, 12-13 uur
18 nov : 16 dec: 20 jan:
“Statistische modellen” “Interim analyses” “”
Sprekers:
Sacha la Bastide, Hans Burgerhof, Václav Fidler afdeling Epidemiologie 2
Voorbeeld: effect van dieet op bloeddruk
Effect van dieet op bloeddruk: analyse proefopzet: 2 onafhankelijke groepen
dieet (A of B)
bloeddruk
groep aantal A n1 B n2
gemiddelde ...
SD ...
...
...
randomizatie confounders
Statistische analyse: t-toets & betrouwbaarheidsinterval voor het verschil tussen gemiddelden toets van Wilcoxon & b.i. voor het verschil tussen medianen
proefopzet: 2 onafhankelijke groepen groep aantal A n1 B n2
gemiddelde ...
SD ...
...
...
3
4
Heb ik een statistisch model nodig?
Bloeddrukdata: statistische model
Statistisch model groep A B
waarnemingen X1,...,Xn Y1,...,Yn
kansverdeling
Statistisch model: Vertaling van proefopzet (en vraagstelling) in termen van kansverdeling van waarnemingen.
Xi ~ kansverdeling: N(µA, σ2) Yi ~ N(µB, σ2) H0: µA = µB
(onafhankelijk)
Xi ~ kansverdeling: FA(u)=P(X
P(YY)>½ of P(X >Y)<½
5
6
1
Stochastisch groter
Voorbeeld 2: leeftijd en hypertensie
1
leeftijd
hypertensie (ja / nee)
proefopzet: observationeel nagenoeg dezelfde situatie als in vorig voorbeeld 2 onafhankelijke groepen ...
Y X
hypertensie aantal gemiddelde ja n1 ... nee n2 ... Statistische analyse: t-toets, Wilcoxon ...
P(X P(Y
SD ... ...
bloeddruk u 7
echter: niet handig voor causale interpretatie beter: P(hypertensie gegeven de leeftijd)
8
Data van 400 personen boven 60 jaar
Leeftijd en hypertensie (2) P(hypertensie) = P(Y=1) = functie van leeftijd P(Y = 1) = π = logit(π )=log(
ea +b⋅leeftijd : logistische regressie e a +b⋅leeftijd + 1
π
1− π
) = a + b ⋅ leeftijd
P(Y = 1) = π = Φ (a + b ⋅ leeftijd ), Φ : verdelingsfunctie van de normale verdeling probit(π )=Φ −1 (π ) = a + b ⋅ leeftijd
andere mogelijkheden: Cauchy verdeling, … 9
10
Opdrachten in R
Resultaten opvragen
> plot(dd$age,dd$sbp, xlab="age", ylab="SBP", type = "n“, xlim=c(60,90), xaxp=c(60,90,3), ylim=c(100,250), yaxp=c(100,250,3) ) > points(dd$age,dd$sbp,cex=0.5, col="blue")
> lm1 <- lm( sbp ~ I(age-65), data=dd ) > summary(lm1) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 159.8290 1.2706 125.792 < 2e-16 *** I(age - 65) 0.6348 0.1938 3.275 0.00115 ** --Residual standard error: 23.76 on 398 degrees of freedom Multiple R-squared: 0.02625, Adjusted R-squared: 0.0238 F-statistic: 10.73 on 1 and 398 DF, p-value: 0.001147
“linear model”
> lm1 <- lm( sbp ~ age, data=dd ) > abline(lm1, col="red“ )
R werkt met objecten (getallen, variabelen, data-frames, functies en procedures)
Op welke leeftijd is de gemiddelde bloeddruk 160? Data staan in dataframe dd, met variabelen age, sbp,... aan te roepen met b.v. dd$age
65 + (160-159.829)/0.6348 = 65.27 Boven deze leeftijd is de gem.bloeddruk hoger dan 160. 11
12
2
Hypertensie: dichotoom SBPh = 1 =0
Hypertensie: t-toets > t.test( dd$age[dd$sbph==0], dd$age[dd$sbph==1] )
als SBP>160 anders
-
t-toets
-
logistische regressie probit-analyse kwantiel-regressie
Welch Two Sample t-test t = -2.6039, df = 382.085, p-value = 0.009575 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -2.7921799 -0.3896383 sample estimates: mean of x mean of y 66.62500 68.21591
13
De toets is correct, maar het resultaat is in deze vorm moeilijk te interpreteren.
Logistische regressie
14
Logistische regressie (2)
> lm2 <- glm( sbph ~ I(age-65), data=dd, family=binomial(link = "logit") )
Coefficients: Estimate Std. Error z value Pr(>|z|)
> summary(lm2) Coefficients:
(Intercept) -0.34287
0.10933
-3.136
I(age - 65)
0.01665
2.554
0.04253
0.00171 ** 0.01064 *
Estimate Std. Error z value Pr(>|z|) (Intercept) -0.34287
0.10933
-3.136
I(age - 65)
0.01665
2.554
0.04253
0.00171 **
Odds ratio voor hypertensie gekoppeld aan leeftijdsverschil van 10 jaar is exp(0.4253)=1.53.
0.01064 *
Op welke leeftijd is de kans op hypertensie meer dan 50%? 65 + 0.34287/0.04253 = 73.06 jaar 15
Probit-analyse
Logit en probit - analyse
> lm3 <- glm( sbph ~ I(age-65), data=dd, family=binomial(link = “probit") )
- Modellen verschillen niet veel.
> summary(lm3)
- Coefficienten in logistische regressie zijn te interpreteren in termen van odds-ratios.
Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.21474
0.06788
-3.164
I(age - 65)
0.01034
2.566
0.02654
16
- Coefficienten in probit-analyse zijn (moeilijker) te interpreteren in de termen van de normaalverdeelde lineaire predictor.
0.00156 ** 0.01028 *
- Dichotomiseren betekent informatie-verlies. Resultaten kunnen afhankelijk zijn van de gekozen cut-off.
Op welke leeftijd is de kans op hypertensie meer dan 50%? 65 + 0.21474/0.02654 = 73.09 jaar 17
18
3
Mediane-regressie
Kwantiel-regressie
Normale lineaire regressie minimaliseert de som van de kwadratische afwijkingen:
∑(y − a −b⋅ x ) i
minimum van
i
Mediane regressie minimaliseert de som van absolute afwijkingen:
i
wordt bereikt voor a = τ − kwantiel
0 < tau < 1 Kwantiel-regressie minimaliseert
∑| y − a − b ⋅ x | i
∑ ρτ ( y − a)
waaarbij ρτ ( x) = τ ⋅ x ⋅ I ( x > 0) + (τ − 1) ⋅ x ⋅ I ( x < 0)
2
∑ ρτ ( y
i
−a − b ⋅ x i )
i
tau=0.5 : mediane regressie 2
minimum van
∑ ( yi − a) wordt bereikt voor a = Y
minimum van
∑| y − a | i
wordt bereikt voor a = mediaan 19
Boxplot van SBP per leeftijdsklasse
20
Mediane regressie > library(quantreg) > lm3 <- rq(sbp ~ I(age-65), tau=0.5, data = dd) > summary(lm3)
240 220
Coefficients: Value Std. Error t value Pr(>|t|) (Intercept) 159.00000 1.19705 132.82610 0.00000 I(age - 65) 0.33333 0.23247 1.43388 0.15239
200 180 160
Op welke leeftijd is de mediane bloeddruk 160?
Modeleer (een aantal) kwantielen in plaats van alleen het gemiddelde.
140 120
(60,63]
(66,69]
(72,75]
(78,81]
65 + (160-159)/0.3333 = 68.3 Boven deze leeftijd is de mediane bloeddruk hoger dan 160.
(84,87]
21
Kwantiel-regressie tau=(0.05, 0.1, 0.25, 0.5, 0.75, 0,9, 0.95)
22
Kwantiel-regressie: 49 kwantielen met age, sex, smoking, BMI
250
> plot(summary(rq(sbp ~ I(age-65)+…,tau = 1:49/50,data=dd)))
mean (LSE) fit median (LAE) fit
SBP
200
150
100 60
70
80 age
90 23
24
4
Kwantiel-regressie
Conclusie - Bij het vergelijken van twee onafhankelijke groepen hoeft de t-test niet de meest nuttige analyse zijn. Formuleren van het onderliggende model helpt. - Bij twijfel over de veronderstellingen van de lineaire regressie (normaliteit, constante variantie) kan men naast het transformeren of dichotomiseren ook de kwantiel-regressie toepassen.
25
26
Volgende Help! Statistiek! lezing:
woensdag 16 december 2009, 12-13 uur zaal 16, UMCG Onderwijs Centrum
“Interim analyses”
Presentatie komt te staan in Download area op http://www.EpidemiologyGroningen.nl
27
5
