kvadrater = rep(0, 10) # Skapar vektor med nollor där vi sen sparar varje kvadrat
for (i in 1:10){
kvadrater[i] = i^2
}
kvadrater [1] 1 4 9 16 25 36 49 64 81 100Statistik och dataanalys I, 15 hp
Datorlaboration 6 - Konfidensintervall, hypotesttest och centrala gränsvärdessatsen
Introduktion och viktig förberedande information
I den här datorlabben kommer vi att arbeta med konfidensintervall och hypotestest
Installera paket
Den här labben förutsätter inga installerade paket.
Skapa mapp för labben
Skapa en mapp Lab6 i din kursmapp SDA1. Ladda ner Quarto-filen för denna lab genom att högerklicka här och välj ‘Spara länk’ eller något likande från menyn som dyker upp. Spara filen i din Lab6 mapp. Öppna Quarto-filen i RStudio och fortsätt med denna laboration direkt i Quarto-dokumentet, där du också fyller i svaren på dina laborationsövningar.
Looping as a way of life
I den här uppgiften kommer ni behöva upprepa beräkningar ett stort antal gånger. Då använder vi loopar, framförallt den s k for-loopen. Här är en enkel loop som beräknar kvadraten av alla tal mellan 1 och 10 och sparar resultaten i en vektor kvadrater:
Vi kan göra precis samma loop på lite annorlunda sätt:
kvadrater = rep(0, 10)
min_loop_vektor = 1:10
for (i in min_loop_vektor){
kvadrater[i] = i^2
}Notera att R alltså kan loop:a över vilken vektor som helst. R läser for-loopen ovan som ’upprepa kommandot print(i^2) för alla värden på i tagna från vektorn min_loop_vektor och spara resultatet i den i:te positionen i vektorn kvadrater.
Genom att använda seq() funktionen kan vi t ex loop:a över värdena 0.1, 0.2,…, 0.9, 1 så här:
min_loop_vektor = seq(0.1, 1, by = 0.1)
for (i in min_loop_vektor){
print(i^2)
}[1] 0.01
[1] 0.04
[1] 0.09
[1] 0.16
[1] 0.25
[1] 0.36
[1] 0.49
[1] 0.64
[1] 0.81
[1] 1Men hur gör vi om vi vill spara dessa kvadrater i en vektor kvadrater som vi gjorde ovan? Vi kan inte längre använda loop-variabeln i för tala om för R att placera kvadraten i^2 på den i:te positionen i kvadraten. Variabeln i räknar ju inte längre från 1 till 10! Den räknar ju 0.1, 0.2,…, 0.9, 1.
Lösningen är att skapa ytterligare en variabel som håller reda på var i vektorn kvadrater som resultat ska sparas. Jag kallar den variabel count, men den kan heta vad som helst.
kvadrater = rep(0, 10)
min_loop_vektor = seq(0, 1, by = 0.1)
count = 0
for (i in min_loop_vektor){
count = count + 1
kvadrater[count] = i^2
}
kvadrater [1] 0.00 0.01 0.04 0.09 0.16 0.25 0.36 0.49 0.64 0.81 1.00Notera hur jag först satte count till 0 innan loopen och hur jag sen inne i loopen ökade count med ett i varje upprepning.
Ok, nu testar vi den här kunskapen på ett lite, lite mer avancerat problem. Säg att vi vill simulera från samplingfördelningen för stickprovsandelen \(\hat p\) med data från modellen
\[ X_1,X_2,\ldots,X_n \overset{iid}{\sim}\mathrm{Bernoulli}(p) \]
där p=0.3 och stickprovsstorleken n=100. Vi gör detta genom att simulera 1000 olika stickprov från modellen. För varje stickprov med n=100 observationer beräknar vi andelen \(\hat p\) och sparar denna andel i en lång vektor. Sen plottar vi fördelningen av dessa 1000 st värden på \(\hat p\) med ett histogram. Here we go!
p = 0.3 # andelen i populationen/modellen
nsim = 1000 # vi vill simulera 1000 olika stickprov
n = 100 # varje stickprov består av 100 observationer.
phats = rep(0, nsim) # skapar vektor där andelen (phat) i varje stickprov sparas.
count = 0 # den här variabeln kommer hålla koll på stickprovets nummer
for (i in 1:nsim){
count = count + 1
x = rbinom(n, size = 1, prob = p) # size = 1 gör att binomialfördelningen blir en Bernoulli.
phats[count] = mean(x) # mean(x) ger oss andelen phat.
}
hist(phats, 50)
Som en sista övning med loopar, låt oss beräkna sannolikheten \(P(X\leq 0.5)\) för en Bernoulli-variabel med sannolikheten \(p\) där vi varierar \(p\) mellan 0 och 1 i steg av 0.01:
ps = seq(0, 1, by = 0.01)
probability = rep(0, length(ps)) # funktionen length() räknar antalet element in en vektor
count = 0
for (p in ps){ # man får använda annat namn än i som loop-variabel
count = count + 1
probability[count] = pbinom(0.5, size = 1, prob = p) # p ändras i varje upprepning
}
plot(ps, probability, type = "l", xlab = "p", ylab = "P(X<=0.5)")
1. Konfidensintervall och test för en andel
💪 Uppgift 1.1
På föreläsningen om konfidensintervall för en andel såg vi att andelen S-väljare p skattades till 0.371 i SVT/Novus undersökning i Februari 2023 baserat på ett stickprov av n=3539 personer där alltså 1313 person angav att de skulle rösta på S. Vi såg också att ett 95% konfidensintervall för p var (0.355, 0.387).
Månaden innan (SVT/Novus, Januari 2023) skattades andelen S-väljare till 0.347 (samma antal personer tillfrågades och 1228 personer angav att de skulle rösta på S). Använd funktionen prop.test för att beräkna ett 95% konfidensintervall för p i Januari 2023.
Beräkna även detta konfidensintervall med formeln
\[ \hat p \pm z_{\alpha/2}\cdot\sqrt{\frac{pq}{n}} \]
💪 Uppgift 1.2
I senaste valet fick S 30,33% av rösterna. Använda prop.test för att testa på 5% signifikansnivå om andelen S-väljare i i Februari 2023 är förändrat sedan valet.
2. Konfidensintervall och test för ett väntevärde
I Övning 20.65 i SDA-boken analyseras data för mängden socker i frukostflingor (gram per 100 gram flingor) för flingor som riktar sig till barn (child_sugar):
child_sugar = c(40.3,55.0,45.7,43.3,50.3,45.9,53.5,43,44.2,44,47.4,44,33.6,55.1,48.8,50.4,37.8,60.3,46.6)Låt oss analysera child_sugar med modellen
\[X_1,X_2,\ldots,X_n \overset{iid}{\sim}N(\mu,\sigma)\]
💪 Uppgift 2.1
Använd R för att beräkna ett 95%-igt konfidensintervall för \(\mu\), både genom att använda formeln
\[\bar x \pm t_{\alpha/2,n-1}\frac{s}{\sqrt n}\]
och genom att använda funktionen t.test.
💪 Uppgift 2.2
Använd funktionen t.test för att utföra testet:
\[ \begin{align} &H_0: \mu = 50 \\ &H_1:\mu\neq50 \end{align} \]
genom att beräkna p-värdet för testet och förkasta \(H_0\) om p-värdet är mindre än \(\alpha =0.05\). Dvs testa på 5% signifikansnivå. Tolka resultatet från testet.
💪 Uppgift 2.3
Upprepa testet i uppgiften ovan för olika värden på nollhypotesens värde, \(\mu_0\). Prova alla värden på nollhypotesen \(\mu_0\) i vektorn seq(40, 50, by = 0.1), dvs utför testet i uppgiften ovan med \(\mu_0 = 40\), sen med \(\mu_0=40.1\) osv ända till \(\mu_0=50\). Registera om du förkastar \(H_0\) eller inte för varje värde på \(\mu_0\). För vilka värden på \(\mu_0\) förkastar du \(H_0\)? Jämför med konfidensintervallet i Uppgift 2.1.
Tips
Använd en for-loop för att loop:a över alla värden på \(\mu_0\). Använd min kod ovan som beräkna sannolikheten \(P(X\leq 0.5)\) för en Bernoulli-variabel med sannolikheten \(p\) som en mall för denna loop.
Resultatet från t.test är en lista. Så använd $-tecknet för att plocka ut p-värdet från resultat-listan från testet (t ex resultat$p.value om du sparat resultatet från t.test i variabeln resultat).
💪 Uppgift 2.4
Du är egentligen intresserad av att bevisa att den genomsnittliga sockermängden i barnflingor är mer än 40 gram. Dvs, du vill testa den enkelsidiga hypotesen på signifikansnivån \(\alpha =0.05\):
\[ \begin{align} &H_0: \mu = 40 \\ &H_1: \mu > 40 \end{align} \]
använd funktionen t.test för att testa detta. Tips: läs i hjälpen ?t.test.
3. Centrala gränsvärdessatsen
Antag att data kommer från en likforming (Uniform) fördelning mellan 0 och 1, dvs \(a=0\) och \(b=1\) i den likformiga fördelningen. Alltså: \(X_1,X_2,\ldots,X_n\overset{iid}{\sim}\mathrm{Uniform}(0,1).\)
💪 Uppgift 3.1
Simulera 10000 stickprov med vardera \(n=2\) observationer och beräkna medelvärdet i varje stickprov. Rita ett histogram över dessa 10000 medelvärden.
💪 Uppgift 3.2
Upprepa Uppgift 3.1 för \(n=5\), \(n=10\) och \(n=30\) observationer i varje stickprov. Kommentera resultatet.
💪 Uppgift 3.3
Upprepa Uppgift 3.2, men istället för att beräkna medelvärdet i varje stickprov så ska du beräkna statistikan
\[ x_{\mathrm{max}}=\max(x_1,x_2,\ldots,x_n) \]
dvs det högsta värdet i varje stickprov. Varför verkar centrala gränsvärdessatsen inte fungera här?