Vejledende besvarelse af hjemmeopgave i
Basal statistik for lægevidenskabelige forskere, forår 2013
I forbindelse med reagensglasbehandling blev 100 par randomiseret til to forskellige former for hormonstimulation. Herefter blev der udtaget et eller flere
æg samt (i princippet) to ægblærevæsker, en fra hver æggestok.
Kvaliteten af de udtagne æg blev vurderet, og et antal hormoner blev målt i
ægblærevæskerne, heriblandt testosteron (nmol/l) og InhibinB (ng/ml).
For 4 kvinder lykkedes det ikke at opnå nogle æg eller ægblærer, og for visse
af de 96 resterende var det (af forskellige årsager) ikke muligt at foretage alle
hormonbestemmelserne.
På hjemmesiden
http://staff.pubhealth.ku.dk/~lts/basal13_1/hjemmeopgave.html
ligger oplysninger på 96 kvinder, dels af den anvendte stimulationsmetode
(treat) og dels af de 4 hormonkoncentrationer (to hormoner, på såvel højre
som venstre side).
Datasættet indeholder ikke individuelle registreringer af kvaliteten af de udtagne æg, men nedenfor ses en opgørelse over fordelingen af lav resp. høj
kvalitet af disse, sammenholdt med stimulationsmetoden
Treat=1
Treat=2
Total
Ægkvalitet
dårlig god Total
21
30
51
14
31
45
35
62
96
I de første 4 spørgsmål fokuserer vi udelukkende på målingerne på venstre
side
1. Lav en figur til illustration af forskellene på de to behandlingsgrupper
forsåvidt angår testosteron? Udfør på baggrund af denne tegning en
sammenligning af de to grupper på passende skala og kvantificer forskellen med konfidensinterval.
1
Vi starter med et simpelt scatter plot
proc gplot data=a1;
plot testo_left*treat;
run;
som dog er gjort lidt pænere ved at skrive
title ’spm 1’;
proc gplot data=a1;
plot testo_left*treat
/ frame haxis=axis1 vaxis=axis2;
axis1 offset=(8,8) value=(h=3) minor=NONE label=(H=3);
axis2 value=(h=3) minor=NONE label=(A=90 R=0 H=3);
symbol1 v=circle i=none c=blue h=3;
run;
Man bemærker her, at fordelingen i begge grupper er en anelse skæv,
med en hale mod de høje værdier. Da vi er på vej mod en sammenligning af grupperne med et T-test, forsøger vi at transformere med
2
logaritmen for at opnå symmetriske fordelinger (helst normalfordelinger). Her er anvnedt 10-tals logaritmer, men det spiller ingen rolle.
Ud fra denne figur ser det fornuftigt ud at foretage et uparret T-test,
idet
• Observationerne antages uafhængige (forskellige kvinder)
• Der synes at være nogenlunde samme variation i begge grupper
• Fordelingerne ser rimeligt symmetriske/normalfordelte ud
Man kan selvfølgelig supplere med fraktildiagrammer, men det skønnes
ikke at være nødvendigt her, ikke mindst fordi fordelingsantagelsen ikke er voldsomt kritisk, hvis vi blot skal sammenligne middelværdier og
ikke udtale os om enkeltindivider.
Vi udfører så det uparrede T-test:
proc ttest data=a1;
class treat;
var log_testo_left;
run;
og får outputtet
3
The TTEST Procedure
Variable:
log_testo_left
treat
1
2
Diff (1-2)
treat
1
2
Diff (1-2)
Diff (1-2)
N
48
42
Mean
1.4199
1.3649
0.0550
Std Dev
0.1687
0.1681
0.1684
Pooled
Satterthwaite
Mean
1.4199
1.3649
0.0550
0.0550
Method
Method
Pooled
Satterthwaite
Variances
Equal
Unequal
Std Err
0.0244
0.0259
0.0356
Minimum
1.0175
1.0026
95% CL Mean
1.3709
1.4689
1.3125
1.4173
-0.0157
0.1257
-0.0157
0.1257
DF
88
86.505
t Value
1.55
1.55
Maximum
1.8925
1.7956
Std Dev
0.1687
0.1681
0.1684
Pr > |t|
0.1258
0.1258
Equality of Variances
Method
Folded F
Num DF
47
Den DF
41
F Value
1.01
Pr > F
0.9859
Her ser vi
• Der indgår i alt 90 personer i vores sammenligning, idet nogle
kvinder har manglende værdier for testosteron på venstre side.
• Spredningerne i de to grupper er meget ens (0.1687 hhv. 0.1681,
P-værdi for test af identitet af disse er 0.99.
• Middelværdierne i de to grupper er ikke signifikant forskellige
(P=0.13).
• Den estimerede forskel på de to behandlinger er (på logaritmisk
skala) 0.0550, med et konfidensinterval på (-0.0.0157, 0.1257). Tilbagetransformeres dette, fås et ratio-estimat på 100.0550 = 1.135,
altså at treat=1 ligger 13.5% højere end treat=2, med et konfidensinterval på (0.96, 1.34), dvs. svarende til fra 4% under til 24%
over.
2. Lav en figur til at illustrere sammenhængen mellem inhibinB og testosteron, begge målt på venstre side (benyt evt. symboler svarende til
de to forskellige stimulationer). Udfør på denne baggrund en lineær regressionsanalyse (igen på passende skala), med inhibin (venstre side)
som outcome og testosteron (venstre side) som forklarende variabel. Angiv estimater for intercept og hældning med konfidensintervaller samt
4
spredning omkring regressionslinien. Giv også, så vidt muligt, fortolkninger af disse.
Her er det naturligt at starte med et scatterplot, hvor outcome (inhibin
på venstre side) sættes på Y-aksen, medens kovariaten (testosteron på
venstre side) sættes på X-aksen. Nedenfor er tilføjet regressionslinier
som støtte for vurderingen af rimeligheden af forudsætningerne for at
foretage regressionsanalysen.
Denne figur viser, at residualerne fra en regressionsanalyse på disse
utransformerede data ville have to problemer
• Residualerne har ikke samme spredning for alle niveuer af outcome, men har derimod en tendens til at være større (numerisk, dvs.
både store neagtive og store positive værdier) for de høje (forventede) værdier af outcome.
Dette illustreres også ved det trompetagtige udseende af plottet
af residualer mod forventede værdier for denne regressionsanalyse
på utransformeret skala
5
• Residualerne er ikke normalfordelte, men har en hale mod de store
værdier. Dette illustreres bedst ved et fraktildiagram, der klart har
“hængekøjefacon”
I lyset af disse problemer, vil vi transformere vores outcome med
en logaritme, og vi holder os til 10-tals logaritmen. Der er ikke
noget i de ovenstående problemer, der indikerer, at vi skal transformere vores kovariat. I så fald skulle der fremkomme buer på et
plot af residualer mod kovariaten, men da dette fuldstændig vil
ligne ovenstående plot mod de forventede værdier (bortset fra en
6
spejling), kan vi allerede ud fra det ovenstående se, at der næppe
er nogen problemer med lineariteten.
Alligevel kan det rent fortolkningsmæssigt være rimeligt at transformere kovariaten også, da der i begge tilfælde er tale om koncentrationer. Vi vil derfor her vise begge muligheder, og de forskellige
fortolkninger, de fører til.
Først situationen, hvor kovariaten (testosteron, venstre side) ikke
er logaritmeret:
med den tilhørende regressionsanalyse:
proc glm data=a1;
model log_inhibin_left = testo_left / solution clparm;
run;
The GLM Procedure
Number of Observations Read
Number of Observations Used
96
89
The GLM Procedure
7
Dependent Variable: log_inhibin_left
Source
Model
Error
Corrected Total
R-Square
0.232259
DF
1
87
88
Coeff Var
19.72328
Sum of
Squares
2.41021806
7.96707856
10.37729662
Root MSE
0.302615
Mean Square
2.41021806
0.09157562
F Value
26.32
Pr > F
<.0001
log_inhibin_left Mean
1.534302
Standard
Error
0.07974468
0.00271262
Parameter
Intercept
testo_left
Estimate
1.908854686
-0.013916397
Parameter
Intercept
testo_left
95% Confidence Limits
1.750353501 2.067355871
-0.019308018 -0.008524776
t Value
23.94
-5.13
Pr > |t|
<.0001
<.0001
Vi er nu nede på 89 observationer i analysen, idet der yderligere er
røget en ud på grund af manglende værdier for inhibin på venstre side.
Estimatet for interceptet ses at være 1.909 med tilhørende konfidensinterval (1.750, 2.067), men dette er helt meningsløst, da det refererer
til en testosteronværdi på 0.
Estimatet for hældningen er −0.0139, med et konfidensinterval på
(−0.0193, −0.0085) og dette viser tydeligt den negative association mellem de to hormonværdier. Hvis vi omsætter relationen til den oprindelige skala, har vi
inhibin = 10(1.909−0.0139×testosteron) = 81.068 × 0.97testosteron
Fortolkningen af den tilbagetransformerede hældning 10−0.0139 = 0.97
er altså, at hver gang testosteron øges med 1 enhed, så falder inhibin
med 3%. Hvis testosteron øges med 10 enheder, bliver den tilsvarende
faktor 10−0.0139×10 = 10−0.139 = 0.73, altså et fald på 27%.
Estimatet for spredningen omkring regressionslinien er 0.3026, og dette
skal bruges til konstruktion af normalområder. Man kan godt tilbagetransformere denne spredning, men så bliver det hele meget sværere.
Det er bedst (lettest) at udregne normalområder på logaritmisk skala
og så tilbagetransformere endepunkterne. Dette skal vi gøre i spørgsmål
3 nedenfor.
8
Hvis vi også logaritmetransformerer kovariaten testosteron:
og dernæst laver regressionsanalysen
proc glm data=a1;
model log_inhibin_left = log_testo_left / solution clparm;
run;
finder vi outputtet
The GLM Procedure
Number of Observations Read
Number of Observations Used
96
89
The GLM Procedure
Dependent Variable: log_inhibin_left
Source
Model
Error
Corrected Total
R-Square
0.229120
Coeff Var
19.76355
DF
1
87
88
Sum of
Squares
2.37764809
7.99964853
10.37729662
Root MSE
0.303233
Mean Square
2.37764809
0.09194998
F Value
25.86
log_inhibin_left Mean
1.534302
9
Pr > F
<.0001
Standard
Error
0.03568669
Parameter
log-inhibin for testo 30
Estimate
1.45545953
Parameter
log-inhibin for testo 30
95% Confidence Limits
1.38452837
1.52639069
Standard
Error
0.26656726
0.18965940
Parameter
Intercept
log_testo_left
Estimate
2.879927992
-0.964433621
Parameter
Intercept
log_testo_left
95% Confidence Limits
2.350096689 3.409759295
-1.341402223 -0.587465020
t Value
40.78
t Value
10.80
-5.09
Pr > |t|
<.0001
Pr > |t|
<.0001
<.0001
Estimatet for interceptet ses at være 2.880 med tilhørende konfidensinterval (2.350, 3.410), men dette er næsten ligeså meningsløst som
ovenfor, idet det refererer til en testosteronværdi på 1 (fordi logaritmen
til 1 er 0).
Estimatet for hældningen er −0.9644, med et konfidensinterval på
(−1.341, −0.587), som igen viser den negative association mellem de to
hormonværdier. Hvis vi omsætter relationen til den oprindelige skala,
har vi denne gang
inhibin = 10(2.880−0.9644×log10(testosteron)) = 758.6 × testosteron−0.9644
altså meget tæt på en omvendt proportionalitet.
I dette tilfælde er det således ikke strengt nødvendigt at tilbagetransformere hældningen. Gør man det alligevel, får man effekten af at 10-doble
testosteronværdien, 10−0.9644 = 0.11, altså at inhibin falder med 89%
eller til 11% af hvad den var før.
Interessant er det, at når begge variable (outcome såvel som kovariat)
logaritmetransformeres, så får man samme hældning uanset hvilken
logaritme, man benytter (blot man selvfølgelig benytter samme logaritme på begge). Det betyder, at man også direkte kan se effekten af en
fordobling af testosteronværdien, nemlig en faktor 2−0.9644 = 0.51, på
inhibin-værdien, altså et fald til ca. halvdelen, netop svarende til den
omvendte proportionalitet.
Estimatet for spredningen omkring regressionslinien er 0.3032, næsten
fuldstændigt det samme som ovenfor (0.3026), så det er svært på denne
10
baggrund at sige, hvilken model, der er bedst.
Vi kan lige med et par residualplots (svarende til dem ovenfor for den
utransformerede outcome variable) sikre os, at vores model ser fornuftig
ud.
Det ses, at forudsætningerne ser fornuftige ud, idet residualerne har
nogenlunde samme spredning uanset niveauet af outcome variablen, og
at de er nogenlunde normalfordelte. Samme konklusion fås fra modelkontrol i tilfælde af, at kovariaten ikke er logaritmetransformeret.
11
3. Hvad er den forventede værdi af inhibinB for en kvinde med en testosteron på 30? Og hvad er normalområdet for sådanne kvinder?
Er det usædvanligt at se en kvinde med en inhibinB-værdi på 100 og en
testosteron-værdi på 30?
I modellen, hvor kovariaten testosteron er utransformeret, finder vi den
ønskede predikterede værdi ved hjælp af en estimate-sætning:
proc glm data=a1;
model log_inhibin_left = testo_left / solution clparm;
estimate ’log-inhibin for testo 30’ intercept 1 testo_left 30;
run;
som giver de ekstra linier
Standard
Error
0.03315106
Parameter
log-inhibin for testo 30
Estimate
1.49136278
Parameter
log-inhibin for testo 30
95% Confidence Limits
1.42547146
1.55725411
t Value
44.99
Pr > |t|
<.0001
Den predikterede inhibin værdi er således 101.4955 = 31.0, og normalområdet på logaritmisk skala er 1.4914 ± 2 × 0.3026 = (0.886, 2.097)
Da dette indeholder 2 = log10 (100) er kombinationen ikke usædvanlig. Vi kan også tilbagetransformere normalområdet til 10(0.886,2.097) =
(7.69, 124.91), som jo indeholder 100.
I modellen hvor testosteron også er logaritmetransformeret, skal vi lige
huske at benytte den transformerede værdi i estimate-sætningen, altså
log10(30)=1.477
proc glm data=a1;
model log_inhibin_left = log_testo_left / solution clparm;
estimate ’log-inhibin for testo 30’ intercept 1 log_testo_left 1.477;
run;
12
hvorved vi finder
Standard
Error
0.03568669
Parameter
log-inhibin for testo 30
Estimate
1.45545953
Parameter
log-inhibin for testo 30
95% Confidence Limits
1.38452837
1.52639069
t Value
40.78
Pr > |t|
<.0001
Den predikterede inhibin værdi er her 101.4555 = 28.5, og normalområdet på logaritmisk skala er 1.4555±2×0.3032 = (0.849, 2.062). Da dette
indeholder 2 = log10 (100) er kombinationen ikke usædvanlig. Vi kan også tilbagetransformere normalområdet til 10(0.849,2.062) = (7.06, 115.32),
som jo indeholder 100.
Hvis man havde gennemført disse beregninger i en model, hvor inhibin
ikke var logaritmetransformeret, ville man have fået et normalområde,
der indeholdt negative værdier. Dette bør altid lede til eftertanke om,
hvad der er galt med modellen!!
Bemærk i øvrigt, at det i dette spørgsmål er vigtigt om normalfordelingen og den konstante varians (spredning) holder, idet vi udtaler os
om enkeltobservationer.
I de to næste spørgsmål fokuserer vi på prediktion af inhibin målinger på venstre side. I spørgsmål 2 ovenfor predikterede vi denne
ud fra testosteron på venstre side, men nu inddrager vi også hormoner
målt på højre side som forklarende variable.
4. Er testosteron værdien på højre side ligeså god til at prediktere som
testosteron på venstre side var (i spørgsmål 2)?
Vi springer over scatterplot, modelkontrol osv., da dette helt følger
overvejelserne fra forrige spørgsmål. Endvidere ser vi kun på situationen, hvor kovariaten også er logaritmetransformeret.
13
proc glm data=a1;
model log_inhibin_left = log_testo_right / solution clparm;
estimate ’log-inhibin for testo 30’ intercept 1 log_testo_right 1.477;
run;
som giver outputtet
The GLM Procedure
Number of Observations Read
Number of Observations Used
96
84
The GLM Procedure
Dependent Variable: log_inhibin_left
Source
DF
Sum of
Squares
Model
Error
Corrected Total
1
82
83
1.06811791
9.23914715
10.30726506
R-Square
0.103628
Coeff Var
21.86179
Root MSE
0.335667
Mean Square
F Value
Pr > F
1.06811791
0.11267253
9.48
0.0028
log_inhibin_left Mean
1.535406
Standard
Error
0.04001834
Parameter
log-inhibin for testo 30
Estimate
1.48574765
Parameter
log-inhibin for testo 30
95% Confidence Limits
1.40613843
1.56535687
Standard
Error
0.32279806
0.22805831
Parameter
Intercept
log_testo_right
Estimate
2.522862806
-0.702176814
Parameter
Intercept
log_testo_right
95% Confidence Limits
1.880714634 3.165010978
-1.155857444 -0.248496185
t Value
37.13
t Value
7.82
-3.08
Pr > |t|
<.0001
Pr > |t|
<.0001
0.0028
Med højre sides testosteronværdi som kovariat bliver spredningen omkring regressionslinien estimeret til 0.3357, altså noget større end før
(hvor vi fik 0.3032). Det lyder umiddelbart meget rimeligt, at venstresiden er bedre til at prediktere venstresiden af inhibin.
Bemærk dog lige, at vores sammenligning strengt taget ikke er helt
fair, idet vi i denne sidste analyse er helt nede på 84 observationer i
14
alt. Det skyldes formentlig, at analysen nu involverer målinger fra såvel
højre som venstre side, med deraf følgende større risiko for manglende
værdier.
Vi kan også prøve at lade begge værdier indgå som kovariater på samme
tid:
proc glm data=a1;
model log_inhibin_left = log_testo_left log_testo_right
/ solution clparm;
run;
og finder så:
The GLM Procedure
Number of Observations Read
Number of Observations Used
96
84
Dependent Variable: log_inhibin_left
Source
DF
Sum of
Squares
Model
Error
Corrected Total
2
81
83
2.39257531
7.91468975
10.30726506
R-Square
0.232125
Coeff Var
20.35875
Root MSE
0.312590
Mean Square
F Value
Pr > F
1.19628765
0.09771222
12.24
<.0001
log_inhibin_left Mean
1.535406
Standard
Error
0.31484332
0.29138361
0.30792042
Parameter
Intercept
log_testo_left
log_testo_right
Estimate
2.867504827
-1.072777562
0.118679793
Parameter
Intercept
log_testo_left
log_testo_right
95% Confidence Limits
2.241065388 3.493944267
-1.652539548 -0.493015575
-0.493985254 0.731344841
t Value
9.11
-3.68
0.39
Pr > |t|
<.0001
0.0004
0.7009
altså en virkelig klar indikation af, hvilken prediktor, der er bedst, da
venstreside målingerne er en stærkt signifiknat prediktor, medens højreside målingerne ikke tilføjer noget signifikant til prediktionsevnen,
når man allerede kender venstreside målingerne.
15
5. Hvis vi også udnytter information om inhibin niveauet på højre side
som forklarende variabel, har vi da nogen glæde af at kende de to testosteron målinger?
Ovenfor så vi, at venstre sides testosteron målinger predikterede bedre
end højre sides, men det virker endnu mere oplagt, at inhibin målingerne på modsatte side skulle prediktere bedst. Vi afprøver dette ved
at inddrage alle tre mulige kovariater i en regressionsmodel:
proc glm data=a1;
model log_inhibin_left = log_inhibin_right
log_testo_left log_testo_right / solution clparm;
run;
og finder derved outputtet
The GLM Procedure
Number of Observations Read
Number of Observations Used
96
84
Dependent Variable: log_inhibin_left
Source
DF
Sum of
Squares
Model
Error
Corrected Total
3
80
83
4.66166828
5.64559678
10.30726506
R-Square
0.452270
Coeff Var
17.30162
Parameter
Intercept
log_inhibin_right
log_testo_left
log_testo_right
Parameter
Intercept
log_inhibin_right
log_testo_left
log_testo_right
Root MSE
0.265650
Mean Square
F Value
Pr > F
1.55388943
0.07056996
22.02
<.0001
log_inhibin_left Mean
1.535406
Estimate
Standard
Error
t Value
Pr > |t|
1.436744809
0.512808580
-0.812436059
0.316057226
0.36777209
0.09043553
0.25184875
0.26398697
3.91
5.67
-3.23
1.20
0.0002
<.0001
0.0018
0.2347
95% Confidence Limits
0.704855018 2.168634600
0.332836147 0.692781014
-1.313631044 -0.311241074
-0.209293590 0.841408043
16
Som vi kan se, har vi i hvert fald ikke meget glæde af at kende testosteronværdien på højre side, så denne udelader vi, hvorved vi får
The GLM Procedure
Dependent Variable: log_inhibin_left
Source
DF
Sum of
Squares
Model
Error
Corrected Total
2
81
83
4.56051352
5.74675154
10.30726506
R-Square
0.442456
Coeff Var
17.34784
Parameter
Intercept
log_inhibin_right
log_testo_left
Parameter
Intercept
log_inhibin_right
log_testo_left
Root MSE
0.266360
Mean Square
F Value
Pr > F
2.28025676
0.07094755
32.14
<.0001
log_inhibin_left Mean
1.535406
Estimate
Standard
Error
t Value
Pr > |t|
1.615976232
0.498532102
-0.606889805
0.33682201
0.08988544
0.18475535
4.80
5.55
-3.28
<.0001
<.0001
0.0015
95% Confidence Limits
0.945806077 2.286146386
0.319688279 0.677375925
-0.974494992 -0.239284619
Vi ser, at der stadig ligger noget information i at kende testosteron
værdien på venstre side. I denne model er residualspredningen 0.2664,
altså en del lavere end modellen uden inhibin højre side som kovariat.
Bemærk, at man ikke bør sammenligne R2 -værdierne for modeller, der
ikke har det samme antal kovariater, da denne altid vil være mindre,
når der er flere kovariater med.
Til sammenligning får vi en residualspredning på 0.2818, hvis vi udelukkende medtager højresidige inhibin målinger som kovariat, se nedenfor:
proc glm data=a1;
model log_inhibin_left = log_inhibin_right
/ solution clparm;
run;
17
The GLM Procedure
Dependent Variable: log_inhibin_left
Source
DF
Sum of
Squares
Model
Error
Corrected Total
1
82
83
3.79498204
6.51228302
10.30726506
R-Square
0.368185
Coeff Var
18.35424
Root MSE
0.281812
Mean Square
F Value
Pr > F
3.79498204
0.07941809
47.78
<.0001
log_inhibin_left Mean
1.535406
Standard
Error
0.13913589
0.08814862
Parameter
Intercept
log_inhibin_right
Estimate
0.5973875565
0.6093414159
Parameter
Intercept
log_inhibin_right
95% Confidence Limits
0.3206019573 0.8741731557
0.4339857205 0.7846971114
t Value
4.29
6.91
Pr > |t|
<.0001
<.0001
For en ordens skyld skal vi lige overbevise os selv om, at den multiple regressionsmodel overhovedet var rimelig, dvs. vi skal udføre de
sædvanlige modelkontroltegninger. Det drejer sig om residualer plottet
mod hver af kovariaterne (her kun vist for de to vigtigste), residualer
plottet mod forventede værdier af outcome samt fraktildiagram over
residualerne.
I de to næste spørgsmål fokuserer vi på sideforskelle for testosteron
6. Beskriv den naturligt forekommende variation i testosteron mellem en
kvindes to æggestokke, evt opdelt efter behandlingsgruppe, hvis dette
skønnes relevant. Kvantificer ved hjælp af limits of agreement på den
relevante skala.
Når vi skal beskrive variationen mellem hormon niveauerne på højre
og venstre side, er det oplagt at tage udgangspunkt i forskellene på de
to sider og kvantificere størrelsen af disse. Først skal man først sikre
sig, at det er fornuftigt at se på differenser, dvs. om størrelsen af disse
er nogenlunde uafhængigt at hormon niveauet. Dette gør vi ved hjælp
af et Bland-Altman plot, som blot er et scatter plot af differenserne
mellem højre og venstre, plottet mod gennemsnittet af de selvsamme
værdier:
18
data a1;
set a1;
dif_testo=testo_right-testo_left;
snit_testo=(testo_left+testo_right)/2;
dif_log_testo=log_testo_right-log_testo_left;
snit_log_testo=(log_testo_left+log_testo_right)/2;
run;
proc gplot data=a1;
plot dif_testo*snit_testo=treat
/ vref=0 lv=33 nolegend frame haxis=axis1 vaxis=axis2;
axis1 offset=(8,8) value=(h=3) minor=NONE label=(H=3);
axis2 value=(h=3) minor=NONE label=(A=90 R=0 H=3);
symbol1 v=circle i=none c=blue h=3 w=3 l=1 r=1;
symbol2 v=star i=none c=red h=3 w=3 l=1 r=2;
run;
19
som giver plottet
Her ses en klar trompetfacon, svarende til, at differenserne klart bliver (numerisk) større ved de højere koncentrationer. Forskellen mellem
højre og venstre er derfor bedre udtrykt i relative termer, altså som
procentvise forskelle, så også her skal vi se på logaritmetransformerede
data:
Her ser det pænere ud, pånær en enkelt outlier, så på denne skala vil
20
vi kvantificere differenserne, enten for hver behandling for sig (det gør
vi ikke her), eller fælles for begge behandlinger. Bemærk, at forudsætningerne om de ens spredninger (Bland-Altman plottet) og normalfordelingsantagelsen (fraktildiagrammet nedenfor) her er meget vigtigere
end i de tidligere spørgsmål (med undtagelse af spørgsmål 3), idet vi
her vil udtale os om enkeltindivider og ikke blot gruppegennemsnit.
For at kvantificere differenserne, udregner vi de relevante summary statistics:
proc means N Nmiss mean std data=a1;
var dif_log_testo;
run;
og får outputtet
The MEANS Procedure
Analysis Variable : dif_log_testo
N
N
Miss
Mean
Std Dev
-----------------------------------------85
11
0.0099032
0.1232715
------------------------------------------
21
Ud fra gennemsnit og spredning af differenserne på logaritmsik skala,
udregner vi limits of agreement:
0.0099032 ± 2 × 0.1232715 = (−0.2366398, 0.2564462),
som tilbagetransformeret giver estimatet 100.0099032 = 1.02 for ratio
mellem siderne, med tilhørende limits of agreement 10(−0.2366398,0.2564462) =
(0.58, 1.80)
Vi kan således konkludere, at 95% af forskellene højre/venstre ligger
mellem ratioerne 0.58 og 1.80, svarende til at højre æggestok ligger
mellem 42% under til 80% over venstre æggestok. Da der næsten ikke
er nogen systematisk forskel (højre æggestok ligger i gennemsnit kun
2% højere end venstre æggestok), ville vi finde ganske det tilsvarende, hvis vi var startet med at se på forskelle ’den anden vej’, altså
venstre-højre. Dette virker måske sært pga de tilsyneladende asymmetriske grænser, men man skal tænke på, at hvis højre side ligger 80%
over venstre side (svarende til f.eks. 20 og 36), så ligger venstre side
16/36 × 100 = 44.4% under højre side.
7. Er der systematisk forskel på testosteronværdierne på højre og venstre
side? Kvantificer på passende skala, med konfidensinterval, og giv en
fortolkning.
Dette spørgsmål besvares ved at udføre et T-test. Dette kan gøres ved
hjælp af proc means, og hvis man bruger option clm, får man også
konfidensintervallet med.
proc means N Nmiss mean stderr clm probt;
var dif_log_testo;
run;
med outputtet
The MEANS Procedure
Analysis Variable : dif_log_testo
22
N
Lower 95%
Upper 95%
N Miss
Mean
Std Error
CL for Mean
CL for Mean Pr > |t|
-------------------------------------------------------------------------85
11
0.0099032
0.0133707
-0.0166858
0.0364923
0.4610
--------------------------------------------------------------------------
Vi finder ingen signifikant forskel (P=0.46).
For at se, om vi burde behandle de to grupper hver for sig, kunne vi se
på et scatter plot af de to sæt differenser:
eller måske endda et uparret T-test af de to (idet forudsætningerne for
et sådant ses at være rimelige udfra ovenstående plot):
proc ttest data=a1;
class treat;
var dif_log_testo;
run;
The TTEST Procedure
Variable:
dif_log_testo
23
treat
1
2
Diff (1-2)
N
45
40
Mean
0.0286
-0.0111
0.0397
treat
1
2
Diff (1-2)
Diff (1-2)
Method
Pooled
Satterthwaite
Method
Pooled
Satterthwaite
Variances
Equal
Unequal
Std Dev
0.1086
0.1363
0.1224
Std Err
0.0162
0.0215
0.0266
Mean
0.0286
-0.0111
0.0397
0.0397
DF
83
74.422
Minimum
-0.2076
-0.5270
95% CL Mean
-0.00405
0.0612
-0.0547
0.0325
-0.0132
0.0926
-0.0140
0.0934
t Value
1.49
1.47
Maximum
0.3113
0.2864
Std Dev
0.1086
0.1363
0.1224
Pr > |t|
0.1396
0.1452
Equality of Variances
Method
Folded F
Num DF
39
Den DF
44
F Value
1.58
Pr > F
0.1444
som viser, at det ikke er særligt påkrævet at opdele (P = 0.14 for identitet af middelværdier for de to sæt differenser, P = 0.14 for identitet
af de tilsvarende varianser).
Testet for “ingen forskel” mellem højre og venstre kan naturligvis også
udføres som et T-test, men nu er der tale om et parret T-test. Dette kan
udføres enten som et one-sample T-test på differenserne (på logaritmisk
skala), eller som et parret T-test på de to sider (igen logaritmetransformeret):
proc ttest;
var dif_log_testo;
run;
proc ttest data=a1;
paired log_testo_left*log_testo_right;
run;
og som output fås nedenstående
The TTEST Procedure
Variable:
dif_log_testo
24
N
85
Mean
0.00990
Mean
0.00990
DF
84
Std Dev
0.1233
95% CL Mean
-0.0167
0.0365
t Value
0.74
Std Err
0.0134
Minimum
-0.5270
Std Dev
0.1233
Maximum
0.3113
95% CL Std Dev
0.1071
0.1452
Pr > |t|
0.4610
The TTEST Procedure
Difference:
N
85
log_testo_left - log_testo_right
Mean
-0.00990
Mean
-0.00990
DF
84
Std Dev
0.1233
95% CL Mean
-0.0365
0.0167
t Value
-0.74
Std Err
0.0134
Std Dev
0.1233
Minimum
-0.3113
Maximum
0.5270
95% CL Std Dev
0.1071
0.1452
Pr > |t|
0.4610
Vi finder (selvfølgelig ganske som ovenfor) ingen signifikant forskel
(P=0.46).
Den estimerede forskel (højre minus venstre) er (på logaritmisk skala)
0.0099, med et konfidensinterval på
(-0.1267, 0.0365). Tilbagetransformeres dette, fås som før ratio-estimatet
på 1.02, nu med 95% konfidensinterval på (0.96, 1.09). Det betyder,
at værdierne på højre side i gennemsnit ligger 2% over de tilsvarende
i venstre side, men usikkerheden på dette estimat viser, at det ligeså
godt kunne dreje sig om 9% over, eller 4% under.
I de to sidste spørgsmål ser vi på ægkvaliteten, som givet i to gange to
tabellen på forsiden.
8. Giv et estimat for sandsynligheden for høj ægkvalitet, for hver af de to
stimulationsmetoder hver for sig. Ser de to estimater forskellige ud?
Først skal vi indlæse disse ekstra data i en to gange to tabel, og udregne de relevante procenter. Vi sætter tabellen op, så behandlingerne
(treat) er rækkerne, og udfaldene (hhv. god(1) og dårlig(0) ægkvalitet) er kolonnerne. De relevante procenter er således rækkeprocenter,
og vi undertrykker derfor de øvrige:
25
data a1;
input treat kvalitet antal;
cards;
1 0 21
1 1 30
2 0 14
2 1 31
;
run;
title ’spm 8’;
proc freq data=a1;
tables treat*kvalitet / nopercent nocol;
weight antal;
run;
Dette giver os output i form af tabellen
The FREQ Procedure
Table of treat by kvalitet
treat
kvalitet
Frequency|
Row Pct |
0|
1|
---------+--------+--------+
1 |
21 |
30 |
| 41.18 | 58.82 |
---------+--------+--------+
2 |
14 |
31 |
| 31.11 | 68.89 |
---------+--------+--------+
Total
35
61
Total
51
45
96
hvoraf vi aflæser, at treatment 1 har en estimeret sandsynlighed for
god ægkvalitet på 58.82%, hvorimod treatment 2 ligger helt oppe på
68.89%, altså noget højere.
9. Lav et test for om kvaliteten af de udtagne æg afhænger af stimulationsmetoden. Kvantificere forskellen, dels som en forskel på de to sandsynligheder fundet ovenfor og dels som en odds ratio. Husk i begge tilfælde
26
at supplere med et konfidensinterval. Kan vi udelukke, at der er dobbelt
så stor sandsynlighed (eller odds) for at have god ægkvalitet i den ene
gruppe i forhold til den anden?
Vi udvider vores analyse af to gange to tabellen med et par ekstra
options:
• expected: De forventede counts i hver celle, under hypotesen om
“ingen forskel på behandlingerne”. dem skal vi bruge til at vurdere
om Chi-i-anden testet er tilladeligt at udføre
• chisq: Chi-i-anden testet skal udskrives
• riskdiff: Estimater for forskellen mellem sandsynlighederne for
god ægkvalitet skal udskrives
• relrisk: Estimat for den relative risiko skal udskrives
title ’spm 8’;
proc freq data=a1;
tables treat*kvalitet / nopercent nocol expected chisq riskdiff relrisk;
weight antal;
run;
som giver outputtet
The FREQ Procedure
Table of treat by kvalitet
treat
kvalitet
Frequency|
Expected |
Row Pct |
0|
1|
---------+--------+--------+
1 |
21 |
30 |
| 18.594 | 32.406 |
| 41.18 | 58.82 |
---------+--------+--------+
2 |
14 |
31 |
| 16.406 | 28.594 |
| 31.11 | 68.89 |
---------+--------+--------+
Total
35
61
Total
51
45
96
27
Statistics for Table of treat by kvalitet
Statistic
DF
Value
Prob
-----------------------------------------------------Chi-Square
1
1.0455
0.3066
Likelihood Ratio Chi-Square
1
1.0506
0.3054
Continuity Adj. Chi-Square
1
0.6561
0.4179
Mantel-Haenszel Chi-Square
1
1.0346
0.3091
Phi Coefficient
0.1044
Contingency Coefficient
0.1038
Cramer’s V
0.1044
Fisher’s Exact Test
---------------------------------Cell (1,1) Frequency (F)
21
Left-sided Pr <= F
0.8918
Right-sided Pr >= F
0.2092
Table Probability (P)
Two-sided Pr <= P
0.1010
0.3960
Column 1 Risk Estimates
(Asymptotic) 95%
(Exact) 95%
Risk
ASE
Confidence Limits
Confidence Limits
------------------------------------------------------------------------Row 1
0.4118
0.0689
0.2767
0.5468
0.2758
0.5583
Row 2
0.3111
0.0690
0.1758
0.4464
0.1817
0.4665
Total
0.3646
0.0491
0.2683
0.4609
0.2687
0.4691
Difference
0.1007
0.0975
-0.0905
0.2918
Difference is (Row 1 - Row 2)
The FREQ Procedure
Statistics for Table of treat by kvalitet
Column 2 Risk Estimates
(Asymptotic) 95%
(Exact) 95%
Risk
ASE
Confidence Limits
Confidence Limits
------------------------------------------------------------------------Row 1
0.5882
0.0689
0.4532
0.7233
0.4417
0.7242
Row 2
0.6889
0.0690
0.5536
0.8242
0.5335
0.8183
Total
0.6354
0.0491
0.5391
0.7317
0.5309
0.7313
Difference
-0.1007
0.0975
-0.2918
0.0905
Difference is (Row 1 - Row 2)
Estimates of the Relative Risk (Row1/Row2)
Type of Study
Value
95% Confidence Limits
----------------------------------------------------------------Case-Control (Odds Ratio)
1.5500
0.6677
3.5982
Cohort (Col1 Risk)
1.3235
0.7677
2.2817
Cohort (Col2 Risk)
0.8539
0.6312
1.1551
28
Sample Size = 96
Vi ser af såvel χ2 -test og Fishers eksakte test, at der ikke er signifikant
forskel på ægkvaliteten for de to behandlinger (P=0.31 hhv. 0.40). De
forventede antal er alle et pænt stykke over de krævede 5 (minimum
16.4), så det er ikke nødvendigt at lave Fishers eksakte test, da chi-ianden approksimationen er god nok.
Den estimerede forskel på de to sandsynligheder for god ægkvalitet (Column 2 risk) er 0.1007 (=0.6889-0.5882), altså 10 procentpoint, med
konfidensgrænser (-0.09,0.29). Der kan altså med rimelighed være op
til 29% større sandsynlighed for god ægkvalitet for treatment 2.
Odds ratio estimeres til 1.550 med konfidensinterval (0.667, 3.598), og
dette indeholder et 2-tal. Det samme er tilfældet for den relative risiko for dårlig ægkvalitet (Col1 Risk), men ikke for god ægkvalitet
(Col2 Risk). Man kan derfor ikke sige, at sandsynligheden for god ægkvalitet kan være dobbelt så stor i den ene gruppe i forhold til den
anden, men man kan sige, at den ene gruppe (gruppe 2) kan have dobbelt så store odds for god ægkvalitet.
29