1. No category

Matematikcentrum
Matematisk statistik
Lunds universitet
MASB11 HT14, lp2
Laboration 4
Statistiska test 2015-01-09
Del I: Styrkefunktion
Del II: Standardtest
Syftet med laborationen är att ni ska bekanta er med olika typer av hypotestest och dess
egenskaper:
 Ni ska träna på de grundläggande begreppen inom hypotesprövning (t.ex.
signifikansnivå och styrka) samt vilka slutsatser man kan dra från analysen (del I av
handledningen).
 Ni skall bekanta er med lite av de funktioner som finns i R vad det gäller olika
grundläggande statistiska test, främst parametriska men i viss mån även ickeparametriska test (del II av handledningen). Till varje parametriskt test ges också
motsvarande icke-parametriska alternativ.
Del I: Grundläggande begrepp vid hypotesprövning
Denna del av laborationen fokuserar på de grundläggande egenskaperna hos ett hypotestest,
felrisker, signifikansnivå samt styrka. Det finns också tre videoinspelningar på kursens
hemsida som presenterar begreppen.
Exempel: Muntorrhet
Läkemedel kan ge en nedsatt salivkörtelfunktion, vilket är en riskfaktor för karies och andra
sjukdomar i munhålan. På 7 slumpmässigt utvalda patienter som alla fick samma medicin
mätte man under 5 minuter den så kallade tuggstimulerade saliven. Normal mängd saliv under
dessa förhållanden är 1 ml/min och muntorrhet anses föreligga när mängden saliv understiger
0.7 ml/min. Som modell antog man att salivmängden är normalfördelad med väntevärde µ
och standardavvikelse σ, där σ anses vara 0.5 ml/min. Intressanta frågeställningar är t.ex.:
 Stöder data vår misstanke att medicinen sänker salivproduktionen?
 Om medicinen ger upphov till en genomsnittlig salivproduktion på 0.8 ml/min, hur
troligt är det att vi kommer att missa den nedsatta salivproduktionen med vårt test?
 Hur många patienter ska vi mäta på om vi vill att testet ska upptäcka en nedsatt
salivproduktion på 0.7 ml/min med sannolikheten 0.95?
På kursens hemsida hittar ni data i filen Saliv.RData. Kopiera också ner rutinerna Hypotes.R
och Styrkefkn.R från kursens hemsida. Kortfattade svar till frågorna som ställs i uppgifterna
finns i slutet på Del I av handledningen.
1. De grundläggande begreppen med rutinen Hypotes.R
a) Först vill man undersöka om data från de 7 patienterna stöder vår misstanke att
medicinen sänker salivproduktionen Ställ upp lämpliga hypoteser.
b) Beräkna medelvärdet av mätningarna.
c) Använd rutinen Hypotes.R för att illustrera testets kritiska område då testet utförs på
signifikansnivå α = 0.05. Det aktuella kommandot är hypotes1(σ, n, µ0, α, H1-riktn),
så om hypoteserna är H0: µ = 1; H1: µ < 1 blir kommandot hypotes1(0.5, 7, 1, 0.05,
’<’).
d) Rutinen markerar det kritiska området och anger ett värde k som är gränsen till
området. Hur har k beräknats?
e) Använd ditt beräknade medelvärde för att utföra testet. Vad är din slutsats om H0?
f) Vad är din konkreta tolkning av signifikansnivån α = 0.05 i detta exempel?
g) Undersök hur det kritiska området (gränsen k) ändras då du ändrar signifikansnivån
till α = 0.01. Vad blir din slutsats nu?
2
2. Testets styrka och styrkefunktion med rutinerna Hypotes.R och
Styrkefkn.R
a) Antag nu att genomsnittlig salivutsöndring i riskgruppen i själva verket är 0.8. Då är
förstås H0: µ = 1 falsk och vi vill att vårt test ska upptäcka detta och förkasta denna
hypotes till förmån för hypotesen H1: µ < 1. Sannolikheten att testet verkligen klarar
av detta kallas för testets styrka i punkten 0.8. Använd rutinen Hypotes.R för att
illustrera testets styrka i punkten 0.8. Kommandot är nu
hypotes1(σ, n, µ0, α, H1-riktn, sant µ), så i detta fall skriver du
hypotes1(0.5, 7, 1, 0.05, ’<’, 0.8).
Rutinen ger dig ytterligare en figur som, förutom signifikansnivån α (felrisk av typ I),
även visar β (felrisk av typ II). Vad är den konkreta tolkningen av β i detta exempel?
Hur hänger β ihop med testets styrka?
b) Mer generellt, testets styrka i punkten µ, är P(H0 förkastas då µ är verklig
genomsnittlig salivutsöndring i riskgruppen). Observera att styrkan beror på värdet µ.
I detta exempel gäller att ju mindre µ är i förhållande till µ0 = 1 desto större är chansen
att testet ska upptäcka att H0 inte gäller. Därför är det intressant att studera styrkan
som en funktion av µ, denna funktion betecknas ofta S(µ), (eller h(µ)).
Rutinen Styrkefkn.R ritar upp styrkefunktionen, kommandot är styrka1(σ, n, µ0, α, H1riktn, sant µ) om du vill rita upp funktionen och markera ett speciellt µ-värde. Använd
alltså kommandot styrka1(0.5, 7, 1,0.05, ’<’, 0.8), vilket ger dig den tidigare figuren
plus styrkan som en funktion av µ.
Uppskatta utifrån styrkefunktionen hur stor sannolikheten är att vi med vårt test
kommer att upptäcka att en grupp som bör klassas som muntorra (µ= 0.7) har en sänkt
salivproduktion. Hur många patienter bör vi mäta på om vi med sannolikheten 0.95
verkligen ska upptäcka att muntorra har en sänkt salivproduktion?
Svar:
1a) H0: µ = 1; H1: µ < 1
1d) k=µ0 - λα σ/√n = 1-1.6445 ∙0.05/√7=0.689
1e) Eftersom medelvärdet < 0.689 förkastas H0 på nivå 0.05
1f) Det är 5% risk att vi påstår att riskgruppen har en sänkt salivproduktion när den i själva
verket är normal
1g) H0 kan ej förkastas på nivå 0.01
2a) β=P(ej förkasta H0 då verklig genomsnittlig salivproduktion i riskgruppen är 0.8). β=1S(0.8), d.v.s. β är 1- styrkan i punkten 0.8
2b) Styrkan i punkten 0.7 är S(0.7) vilket enligt figuren kan uppskattas till 0.48. Det krävs
n=30 patienter för att styrkan ska vara 0.95 i punkten 0.7.
3
Del II: Statistiska standardtest i R
Kortfattade svar till frågorna som ställs i uppgifterna finns i slutet på handledningen.
På kursens hemsida hittar ni de fem datamaterial som används i denna del av handledningen.
Datamaterialen som används är Albumin, Median, Klorofyll, Dammar och Gråsparvar.
1
Test av medelvärde/median från en population.
Använd datamaterialet Albumin.RData.
a) t-test. En blandning av blodserum innehåller exakt 42 g albumin per liter. Två
laboratorier (A och B) får göra sex bestämningar var av koncentrationen. Vi vill
undersöka om det finns någon systematisk avvikelse från det sanna värdet (42 g/l)
(tvåsidig mothypotes). I R gör man två t-test med nedanstående kommandon i
kommandofönstret (glöm inte att spara det som, t.ex. lab3.R) och kör dem.
t.test(AlbuminA, mu=42)
t.test(AlbuminB, mu=42)
Tolka utskriften. Hur stora är p-värdena och vad blir slutsatserna? Vad blir
konfidensintervallen för medelkoncentrationerna?
b) Test av median. Tillämpning av Binomialfördelningen. Använd datamaterialet
Median.RData. I materialet finns total längd för ett stickprov om 16 gråsparvar. Vi vill
testa om populationens medianlängd är större än 160 mm. Detta kan göras genom att
undersöka hur många gråsparvar i stickprovet som är längre än 160 mm. Under H0 är
dessa antal Bin(n, 0.5). Detta innebär att vi kan använda ett test för p i en
binomialfördelning. Vi behöver då n (=16) och beräkna x = antal som är längre än
160 mm och testa H0: p=0.5:
antal <- length(Median)
stora <- sum(Median > 160)
binom.test(x=stora, n=antal, p=0.5)
Lägg märke till att det p-värde vi får gäller ett tvåsidigt test. Om vi vill ha ett ensidigt test
får vi använda ett av kommandona
binom.test(x=stora, n=antal, p=0.5, alternative=”less”)
binom.test(x=stora, n=antal, p=0.5, alternative=”greater”)
Om vi vill testa om medianen är större än 160, vilket alternativ ska vi då använda?
Vad blir resultatet om vi vill ha ett ensidigt test? Hur stort är p-värdet och vad blir
slutsatsen?
c) Rang-test. Alternativ variant av mediantestet är att använda ranger istället. När vi bara
undersöker hur stor andel av värdena som ligger över resp. under den önskade medianen
förlorar vi informationen om hur långt från medianen de ligger. Om de värden som ligger
över 160 dessutom ligger längre från 160 än vad de värden som ligger under 160 gör tyder
det ju på att medianen borde vara större än 160. Detta utnyttjar man i ett Wilcoxon signed
rank-test:
4
wilcox.test(Median, mu=160)
wilcox.test(Median, mu=160, alternative=”less”)
wilcox.test(Median, mu=160, alternative=”greater”)
Vad blir slutsatsen nu? Om vi vill testa om medianen är större än 160, vilket ensidigt
alternativ ska vi välja? (Ignorera varningarna om cannot compute exact test with ties.)
2
Test vid jämförelse av två populationer.
Använd datamaterialet Klorofyll.RData.
a) t-test (oberoende grupper). Alger fick växa under ljusa respektive mörka
förhållanden. Undersök med ett t-test om det finns skillnader i förväntad klorofyllhalt
mellan de två grupperna genom att använda t.test igen men nu genom att lägga till
gruppvariabeln (jämför Lab1 där vi delade upp överlevnaden på hannar och honor
med ~):
t.test(Klorofyll$Klorofyll ~ Klorofyll$Grupp)
Tolka utskriften. Om vi inte säger något annat förutsätter t.test att varianserna i de två
grupperna är olika och kompenserar för det. Om vi vet (eller antar) att varianserna är
lika kan vi utnyttja det:
t.test(Klorofyll$Klorofyll ~ Klorofyll$Grupp, var.equal=TRUE)
Blev det någon skillnad? För säkerhets skull är det bäst att verkligen testa om
varianserna är olika:
var.test(Klorofyll$Klorofyll ~ Klorofyll$Grupp)
Är varianserna olika? Vilket av de två t-testen ska vi använda? Slutsatser?
Vi kan få en grafisk beskrivning av skillnaderna med en boxplot:
plot(Klorofyll$Klorofyll ~ Klorofyll$Grupp)
Verkar det stämma med slutsatsen ovan?
b) Rangsummetest. Gör nu om testet i 2(a) som ett rangsummetest genom
wilcox.test(Klorofyll$Klorofyll ~ Klorofyll$Grupp)
Hur stort är p-värdet och vad blir slutsatsen? Jämför resultatet med t-testet som ni gjorde i
(a)-uppgiften.
3
Test vid matchade data.
Använd datamaterialet Dammar.RData.
a) t-test. Jämför om det finns någon skillnad i kvävebelastning mellan vår och sommar. Vi
måste nu tala om för t.test att vi har två variabler som är parade:
t.test(Dammar$N_belast_V, Dammar$N_belast_S, paired=TRUE)
5
Finns det några skillnader mellan vår och sommar?
b) Teckentest. Gör nu motsvarande analys som ett icke-parametriskt test med ett teckentest
som räknar hur ofta vår-värdet överstiger motsvarade sommar-värde:
antal <- nrow(Dammar)
stora <- sum(Dammar$N_belast_V > Dammar$N_belast_S)
binom.test(x=stora, n=antal, p=0.5)
Finns det några skillnader? Jämför resultatet med t-testet som ni gjorde i (a)-uppgiften.
c) Rangsummetest. Utnyttja nu storleken på värdena med ett rangsummetest istället:
wilcox.test(Dammar$N_belast_V, Dammar$N_belast_S, paired=TRUE)
Vad blir slutsatsen nu?
4
Test av kategoridata.
Använd datamaterialet Gråsparvar.RData.
a) χ2-test. Vi vill nu testa om vi har en jämn könsfördelning i materialet. Om man gör en
enkel tabell på variabeln kön får man följande resultat:
antal <- table(Grasparvar$Kön)
andel <- prop.table(antal)
cbind(antal, andel)
antal andel
Hane 87
0.6397059
Hona 49
0.3602941
Vi kan utföra testet på flera olika sätt. För det första kan vi göra et χ2-test där vi jämför
den observerade fördelningen med den förväntade under hypotesen om jämn
könsfördelning. Detta görs med
chisq.test(antal)
Hur stort är p-värdet och vad blir slutsatsen?
Om man inte säger något annat testar R om andelarna är lika. Vill man t.ex. testa om
det kan vara 60 % hanar och 40 % honor får man göra så här:
chisq.test(antal, p=c(0.6,0.4))
Vad blir slutsatsen nu?
b) Test av ett proportionstal. Vi kan också göra samma test med hjälp av
binomialfördelningen:
antal <- nrow(Grasparvar)
hannar <- sum(Grasparvar$Kön=="Hane")
binom.test(x=hannar, n=antal, p=0.5)
6
Gör testet. Ändrar sig slutsatsen? En fördel här är att vi på köpet fick ett
konfidensintervall för andelen. Men å andra sidan fungerar det bara när vi har två
kategorier. Har vi fler måste vi använda χ2-testet.
c) Test av två proportionstal. Testa nu om andelen som överlever är lika stor bland honor
som bland hanar. Detta görs med med hjälp av ett χ2-test på en korstabell:
tabell <- table(Grasparvar$Kön, Grasparvar$Överlevnad)
prop.table(tabell, margin=1)
chisq.test(tabell)
Hur stort är p-värdet och vad blir slutsatsen?
Svar:
1. a) A: p=0.033; B: p=0.081; Vi kan påvisa en skillnad för A men inte för B.
95% KI – A: (42,06, 42.94); B: (35.69, 42.51)
b) p=0.038; Vi kan påvisa att Md>160
c) p=0.00024; Vi kan påvisa att Md>160.
2. a) Varianserna är inte olika (p=0.48); t-test: p=0.009; Vi kan påvisa en skillnad i klorofyll.
b) p=0.015; Vi kan påvisa en skillnad. Samma slutsats men inte lika signifikant.
3. a) p=0.262 Vi kan inte påvisa någon skillnad
b) p=0.070; Vi kan inte påvisa någon skillnad
c) p=0.195; Vi kan inte påvisa någon skillnad.
4. a) p=0.0011; Vi kan påvisa en avvikelse från jämn könsfördelning; p=0.345: könsfördelningen
skulle kunna vara 60/40.
b) p=0.0014: Vi kan påvisa en avvikelse från jämn könsfördelning
c) p=0.112: Vi kan inte påvisa någon skillnad i andelen som dör
7