Werkcollege statistiek Thema 1
2024-09-30
In dit COO leer je gemiddelden en (95%) betrouwbaarheidsintervallen berekenen, en absorptie data plotten tegenover de concentraties van een antibioticum.
- Voorbereiding
- Betrouwbaarheidsintervallen (stap voor stap)
- 1.1 Data inlezen
- 1.2 Beschrijvende statistiek
- 1.3 Betrouwbaarheidsinterval voor een gemiddelde
- 1.4 Betrouwbaarheidsinterval voor het verschil tussen twee gemiddelden
- OD data inlezen en analyseren
0. Voorbereiding
Zorg er voor dat je de benodigde databestanden hebt gedownload en lokaal hebt opgeslagen:
Het is ook handig om je ‘working directory’ te veranderen in deze map met databestanden. Dit kun je doen via het menu: Session > Set Working Directory > Choose Directory… > selecteer de map waarin de databestanden staan. Let op tijdens het instellen van je ‘working directory’ zie je de bestanden meestal niet staan.
Verder maken we gebruik van functies uit de psych
package. Als je die nog niet hebt geïnstalleerd, doe dat nu met de
functie install.packages(“psych”). Daarna kunnen we hem
inladen:
library(psych)1: Betrouwbaarheidsintervallen in R
In de zelfstudie is gekeken naar de optische densiteit (OD) bij een concentratie van het antibioticum van 250 mg/l. We gebruiken hier hetzelfde voorbeeld en kijken stap-voor-stap hoe we de data kunnen beschrijven, en hoe een 95% betrouwbaarheidsinterval (BHI) in R te laten berekenen.
1.1 Data inlezen
De data zijn opgeslagen in een het bestand OD250.csv. We
lezen de data in en noemen de data frame d1. We bekijken de
eerste regels met head:
d1 <- read.csv("OD250.csv")dim(d1)## [1] 72 2head(d1)## stam OD250
## 1 1 0.091
## 2 2 0.369
## 3 3 0.488
## 4 4 0.563
## 5 5 0.313
## 6 6 0.3591.2 Beschrijvende statistiek
Vraag 1. Maak de side-by-side boxplot uit de zelfstudie van de OD bij 250 mg/l per stam.
Om een gemiddelde OD en standaarddeviatie per stam te krijgen, kunnen
we steeds een selectie maken op de verschillende stammen en vragen om
die kengetallen. Veel makkelijker: we kunnen de
describeBy() functie gebruiken (we krijgen dan veel meer
dan mean en SD):
# Kerngetallen berekenen per stam
# skew = FALSE zorgt ervoor dat kerngetallen die ons niet interesseren niet geprint worden
groepen <- d1$stam
describeBy(d1$OD250, group = groepen, skew = FALSE) ##
## Descriptive statistics by group
## group: 1
## vars n mean sd median min max range se
## X1 1 12 0.12 0.11 0.09 0.07 0.45 0.38 0.03
## --------------------------------------------------------------------------------------
## group: 2
## vars n mean sd median min max range se
## X1 1 12 0.28 0.11 0.24 0.14 0.49 0.36 0.03
## --------------------------------------------------------------------------------------
## group: 3
## vars n mean sd median min max range se
## X1 1 12 0.46 0.18 0.49 0.06 0.73 0.67 0.05
## --------------------------------------------------------------------------------------
## group: 4
## vars n mean sd median min max range se
## X1 1 12 0.61 0.12 0.58 0.42 0.82 0.4 0.03
## --------------------------------------------------------------------------------------
## group: 5
## vars n mean sd median min max range se
## X1 1 12 0.18 0.15 0.11 0.09 0.59 0.5 0.04
## --------------------------------------------------------------------------------------
## group: 6
## vars n mean sd median min max range se
## X1 1 12 0.33 0.1 0.36 0.12 0.46 0.34 0.03Vraag 2. Welke stam heeft de hoogste gemiddelde OD? Welke stam heeft de kleinste standaarddeviatie? En de grootste?
1.3 Betrouwbaarheidsinterval voor een gemiddelde
Laten we nu kijken naar stam 2. We hebben al een gemiddelde OD uitgerekend, wat is het 95% BHI voor de gemiddelde OD bij een concentratie van het antibioticum van 250 mg/l voor stam 2?
Vraag 3. Voor je het BHI uitrekent, kijk nogmaals naar de boxplot. Mogen we voor zo’n steekproef een BHI uitrekenen? Wat nemen we aan wanneer we dat doen?
We gaan het toch doen. Om een 95% BHI (op basis van de t-verdeling)
uit te rekenen, maken we gebruik van de t-test() functie.
De t-toets wordt in Thema 3 behandeld, hier gebruiken we de functie
alleen voor het 95% BHI. We selecteren eerst een subset van de dataframe
door enkel de rijen te nemen waarin informatie voor stam 2 staan:
# Selecteer stam 2 data
stam2 <- d1[d1$stam == 2,]
# Bereken het 95% betrouwbaarheidsinterval
# De rest van de output negeren we
t.test(stam2$OD)##
## One Sample t-test
##
## data: stam2$OD
## t = 8.7514, df = 11, p-value = 0.000002755
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.2077084 0.3472916
## sample estimates:
## mean of x
## 0.2775De gemiddelde OD is 0.277, met 95% BHI (0.208 – 0.347).
Vraag 4. Geef een interpretatie van dit BHI.
1.4 Betrouwbaarheidsinterval voor het verschil tussen twee gemiddelden
Stel dat we de optical density bij concentratie 250 mg/l willen
vergelijken tussen stammen twee en drie. Dan hebben we twee groepen van
metingen die we hier voor het gemak beschouwen als onafhankelijk van
elkaar (zie overpeinzingen in de zelfstudie). Is het verschil dat we in
de boxplot zien tussen de twee stammen een weerspiegeling van een
werkelijk verschil in de populaties (alle mogelijke OD metingen bij 250
mg/l for stammen 2 en 3), of zou dit verschil door toeval
(steekproefvariabiliteit) kunnen zijn onstaan? Om die vraag te
beantwoorden, kunnen we een 95% BHI voor het verschil in gemiddelden
uitrekenen. Daartoe maken we eerst een selectie van alleen die twee
stammen, en daarna gebruiken weer de t-test() functie.
# We selecteren nu rijen waarin informatie staat
# over stam 2 of (het | teken) stam 3
stam23 <- d1[d1$stam==2 | d1$stam==3,]
# We gebruiken hier een 'r formula' om aan te geven wat de groepen zijn
# in de test (stam23$stam). De waardes als input voor de test zijn de OD waarden (stam23$OD).
# Zie ?t.test voor meer info, onder 'formula'.
t.test(stam23$OD~stam23$stam)##
## Welch Two Sample t-test
##
## data: stam23$OD by stam23$stam
## t = -3.0595, df = 18.418, p-value = 0.006621
## alternative hypothesis: true difference in means between group 2 and group 3 is not equal to 0
## 95 percent confidence interval:
## -0.30930175 -0.05769825
## sample estimates:
## mean in group 2 mean in group 3
## 0.2775 0.4610Het verschil in gemiddelden in de steekproeven is 0.2775-0.461 = -0.184, met 95% BHI (-0.309 – -0.058).
Vraag 5a. Geef een interpretatie van dit BHI. Denk je dat er een verschil is in gemiddelden in de populatie?
Vraag 5b. Wat nemen we aan bij het uitrekenen van dit BHI? Voldoen de gegevens aan deze aannames?
2. OD data inlezen en analyseren
In het bestand ODdatavb.csv staan data van een
experiment zoals je zelf gaat uitvoeren: voor zes verschillende
bacteriestammen kijk je bij verschillende concentraties van een
antibioticum naar de groei(remming). De groei wordt uitgedrukt als OD.
Lees de data in en bekijk het:
OD.dat <-read.csv("OD-eigen-data-vb.csv", sep = ";")OD.dat## Concentratie stam1OD1 stam1OD2 stam1OD3 stam1OD4 stam2OD1 stam2OD2 stam2OD3 stam2OD4 stam3OD1 stam3OD2 stam3OD3
## 1 c10 0.5676 0.5721 0.0731 0.7895 0.8288 0.4495 0.2596 0.6829 0.6073 0.9855 0.8426
## 2 c50 0.2304 0.3814 0.0549 0.2330 0.5397 0.4095 0.3193 0.6058 0.5803 0.9423 0.7905
## 3 c100 0.1103 0.3052 0.0526 0.1397 0.3992 0.3626 0.0490 0.5567 0.4768 0.7816 0.5747
## 4 c250 0.0566 0.0673 0.0569 0.0613 0.2452 0.0586 0.0586 0.5376 0.2569 0.2838 0.2650
## 5 c500 0.0499 0.0520 0.0997 0.0583 0.0745 0.0645 0.0541 0.4871 0.0794 0.1949 0.2907
## 6 NegCont 0.1108 0.0489 0.0547 0.0455 0.0563 0.0489 0.0535 0.4524 0.0535 0.0489 0.4198
## 7 PosCont1 0.6521 0.8136 0.5655 0.7447 0.9293 0.8152 0.6181 0.7486 0.7870 0.9263 0.8449
## 8 PosCont2 0.7286 0.6495 0.7362 0.5945 0.6937 0.8123 0.8311 0.7577 0.8337 0.7921 0.8193
## stam3OD4
## 1 0.66600
## 2 0.62800
## 3 0.66410
## 4 0.59335
## 5 0.42135
## 6 0.04440
## 7 0.77470
## 8 0.78970Let op: Dit bestand komt als het goed is volledig overeen met de structuur op de 96-wells platen die jullie gaan maken. De concentraties antibiotica worden in de rijen aangegeven, en in de kolommen staan steeds 4 herhalingen voor 3 verschillende stammen bacterien.Voor de analyse van je eigen data is het belangrijk dat zowel de variabelenamen (in de eerste rij van de .csv file) als de concentraties in de variabele “concentratie” slechts uit 1 woord bestaat.
Vraag 6. Bekijk de OD-waarden in het dataframe. Wellicht vallen je al een paar dingen op? Zijn er bijzonder hoge of bijzonder lage waarden van OD in de tabel? Welke waarden verwacht je ongeveer voor de negatieve controles? En welke waarden voor de positieve controles?
Omdat we per stam, per concentratie, 4 observaties hebben, kunnen we daar gemiddelden van berekenen, en op die manier met grotere precisie bepalen wat het gemiddelde in de populatie zal zijn voor die stam / concentratie combinatie.
De positieve controles hebben we zelfs in twee rijen, dus 8 keer gemeten. Deze 8 metingen zullen we gebruiken om per stam een gemiddelde met 95% BHI te berekenen, zodat we een goede schatting hebben van het populatie-gemiddelde van de OD waarde in die stam als er geen behandeling met antibiotica zou zijn.
De negatieve controles hebben we maar in 1 rij gemeten, maar in principe zijn alle negatieve controles - hoewel per stam gemeten - allemaal schatters voor hetzelfde, namelijk de OD bij een monster zonder bacteriën. Deze kunnen we dus ook middelen, en met een 95% BHI van het gemiddelde hebben we 95% zekerheid dat het BHI het werkelijke populatiegemiddelde van OD zonder infectie bevat.
Als het goed is zullen de (gemiddelden van) de OD waarden van met antibiotica behandelde bacteriele infecties ergens tussen die positieve en negatieve controles in vallen. Maar dit is een experiment, en daar kan natuurlijk vanalles fout gaan.
We gaan aan de slag met de data! Dat kan echter niet zomaar. De data staan nog niet zo geordend dat we meteen berekeningen uit kunnen voeren. We willen graag dat de variabelen in kolommen weergegeven worden in plaats van in rijen, en daarom moeten we deze data file zo herordenen dat de rijen kolommen worden en vice versa. Dat heet transponeren. De onderstaande code zorgt hiervoor. Daarna maken we nog een extra variabele die aangeeft van welke bacteriestam de data zijn.
# Hier transponeren we de data.
# Ook laten we kolom 1 weg, omdat we die later toevoegen als kolomnamen
tmp.dat <- OD.dat[,-1]
tmp.dat <- t(tmp.dat)
# t() maakt een matrix van je data.frame. Dit moeten we dus weer omzetten
tmp.dat <- as.data.frame(tmp.dat)
# We wijzen kolomnamen toe, we gebruiken de waardes die eerst in kolom 1 stonden.
tmp.dat <- setNames(tmp.dat, OD.dat[,1])
# Nu maken we twee extra kolommen, met info over de stam en welke meeting het was.
# Met cbind voegen we deze kolommen toe.
x <- setNames(as.data.frame(c(1,1,1,1,2,2,2,2,3,3,3,3)), "Stam")
y <- setNames(as.data.frame(c(1,2,3,4,1,2,3,4,1,2,3,4)), "Herhaling")
od.dat <- cbind(x, y, tmp.dat)
# Onze nieuwe dataset
od.dat## Stam Herhaling c10 c50 c100 c250 c500 NegCont PosCont1 PosCont2
## stam1OD1 1 1 0.5676 0.2304 0.1103 0.05660 0.04990 0.1108 0.6521 0.7286
## stam1OD2 1 2 0.5721 0.3814 0.3052 0.06730 0.05200 0.0489 0.8136 0.6495
## stam1OD3 1 3 0.0731 0.0549 0.0526 0.05690 0.09970 0.0547 0.5655 0.7362
## stam1OD4 1 4 0.7895 0.2330 0.1397 0.06130 0.05830 0.0455 0.7447 0.5945
## stam2OD1 2 1 0.8288 0.5397 0.3992 0.24520 0.07450 0.0563 0.9293 0.6937
## stam2OD2 2 2 0.4495 0.4095 0.3626 0.05860 0.06450 0.0489 0.8152 0.8123
## stam2OD3 2 3 0.2596 0.3193 0.0490 0.05860 0.05410 0.0535 0.6181 0.8311
## stam2OD4 2 4 0.6829 0.6058 0.5567 0.53760 0.48710 0.4524 0.7486 0.7577
## stam3OD1 3 1 0.6073 0.5803 0.4768 0.25690 0.07940 0.0535 0.7870 0.8337
## stam3OD2 3 2 0.9855 0.9423 0.7816 0.28380 0.19490 0.0489 0.9263 0.7921
## stam3OD3 3 3 0.8426 0.7905 0.5747 0.26500 0.29070 0.4198 0.8449 0.8193
## stam3OD4 3 4 0.6660 0.6280 0.6641 0.59335 0.42135 0.0444 0.7747 0.7897Om de data goed te beschrijven gaan we kijken naar:
- de OD waardes van de negatieve controles
- de OD waardes van de positieve controles
- de OD waardes van de observaties met verschillende concentraties antibiotica
Beantwoord nu op basis van od.dat de volgende vragen
over de OD-waarden van de negatieve controles.
Vraag 7. Beschrijf (getallen, plaatjes) de OD van de negatieve controles. Wat zie je?
Vraag 8. Wat is het 95% BHI voor de gemiddelde OD van de negatieve controles in dit experiment?
Vraag 9. Wat vind je van het gemiddelde OD voor de negatieve controles?
Vraag 10. Wat vind je van de breedte het 95% BHI van het gemiddelde negatieve controle? Hoe zouden we het interval smaller kunnen krijgen?
Het is handig om het gemiddelde en de uiteinden van het BHI van de
negatieve controle op te slaan voor de plot. We halen de onder- en
bovengrens (elementen [1] en [2] van de conf.int) uit het
object dat gemaakt wordt door de t.test() functie:
# Mean
Negmean.nc <- mean(od.dat$NegCont)
# Confidence interval
Output.ttest <- t.test(od.dat$NegCont) # Je ziet dat dit een list is met 10 elementen
# In het element 'conf.int' zit een vector met twee waardes
Negll.nc <- Output.ttest$conf.int[1] # De lower limit (ll)
Negul.nc <- Output.ttest$conf.int[2] # De upper limit (ul)Hoe kun je checken of dit gelukt is?
Dezelfde vragen gaan we nu beantwoorden over de OD-waarden van de positieve controles. Maar daar is eerst weer een bewerking van de data voor nodig. Per stam en herhaling heb je twee positieve controle metingen gedaan. Deze metingen staan nu nog in twee kolommen. Voor je analyse is het makkelijker als je de positieve metingen allemaal in dezelfde kolom hebt. Dit regel je met “reshape”, een functie die van een ‘wide’, naar een ‘long’ format gaat. Je kunt daarvoor de onderstaande code gebruiken. Let op: ook als ervaren R-gebruiken is dit vaak even uitproberen. Als je dit zelf gaat doen is het niet gek als je dit een paar keer moet proberen, tot dat je de indeling hebt die je wilt.
# Eerst selecteren we alleen de gegevens van stam 1
stam1 <- od.dat[od.dat$Stam==1,]
stam1 # wide## Stam Herhaling c10 c50 c100 c250 c500 NegCont PosCont1 PosCont2
## stam1OD1 1 1 0.5676 0.2304 0.1103 0.0566 0.0499 0.1108 0.6521 0.7286
## stam1OD2 1 2 0.5721 0.3814 0.3052 0.0673 0.0520 0.0489 0.8136 0.6495
## stam1OD3 1 3 0.0731 0.0549 0.0526 0.0569 0.0997 0.0547 0.5655 0.7362
## stam1OD4 1 4 0.7895 0.2330 0.1397 0.0613 0.0583 0.0455 0.7447 0.5945# Om naar long te gaan gebruiken we reshape.
# Bekijk ook ?reshape
od.posS1 <- reshape(stam1, # Onze input df
idvar=c("Stam", "Herhaling"), # De id kolommen voor individuele metingen
varying=c("PosCont1", "PosCont2"), # De kolommen met de waardes die we onder elkaar willen
v.names = "OD", # De naam die we geven aan de nieuwe kolom met waardes onder elkaar
direction="long", # We willen naar een 'long' format
times=c("PosCont1", "PosCont2"), # De twee kolommen waar de metingen van de positieve controles in staan
drop=c("c500", "c250", "c100", "c50", "c10","NegCont") # De kolommen die we niet mee willen nemen
)
od.posS1 # long, vergelijk alle waarden: staan die allemaal op de goede plaats?## Stam Herhaling time OD
## 1.1.PosCont1 1 1 PosCont1 0.6521
## 1.2.PosCont1 1 2 PosCont1 0.8136
## 1.3.PosCont1 1 3 PosCont1 0.5655
## 1.4.PosCont1 1 4 PosCont1 0.7447
## 1.1.PosCont2 1 1 PosCont2 0.7286
## 1.2.PosCont2 1 2 PosCont2 0.6495
## 1.3.PosCont2 1 3 PosCont2 0.7362
## 1.4.PosCont2 1 4 PosCont2 0.5945Nu kunnen we onderstaande vragen beantwoorden.
Vraag 11. Beschrijf (getallen, plaatjes) de OD van de positieve controles voor de eerste stam. Wat zie je?
Vraag 12. Wat is het 95% BHI voor de gemiddelde OD van de positieve controles van stam 1 in dit experiment?
Vraag 13. Wat vind je van het gemiddelde en de breedte het 95% BHI van de positieve controles van stam 1?
Ook het gemiddelde en de uiteinden van het BHI van de positieve controle slaan we op voor de plot:
PosS1mean.nc <- mean(od.posS1$OD)
PosS1ll.nc <- t.test(od.posS1$OD)$conf.int[1]
PosS1ul.nc <- t.test(od.posS1$OD)$conf.int[2]Herhaal bovenstaande voor de positieve controles van stammen 2 en 3.
Als laatste gaan we nu per stam een grafiek maken waarin de OD waarden van verschillende concentraties gevisualiseerd worden.
Om grafieken te maken per stam, starten we weer met de datafile
od.dat, en moeten we de data weer omvormen. De gegevens staan nu weer in
het ‘wide format’, met één regel per stam, maar om een grafiek te
krijgen van OD per concentratie moeten we één regel per concentratie per
stam en per herhaling hebben (‘long format’). We gebruiken dus weer de
reshape() functie. Omdat deze functie meestal wordt
gebruikt bij longitudinale data, heet de optie om de verschillende
regels per concentratie te maken ‘times’. Dat is wellicht verwarrend,
‘time’ is in dit geval dus concentratie. Als we kijken naar de
resulterend data frame, dan zien we dat het handig is om de tweede
variabele te hernoemen. Dat doen we met colnames().
od.long <- reshape(od.dat, # Input df
idvar=c("Stam","Herhaling"), # Variabele om de verschillende meetingen uit elkaar te houden
varying=list(3:7), # De kolommen met OD waardes, we gebruiken nu numerieke indexen ipv kolom namen
v.names = "OD", # Naam van de nieuwe kolom met OD waardes
direction="long",
times=c(10, 50, 100, 250, 500), # Let op dat deze volgorde klopt voor je eigen data
drop=c("NegCont","PosCont1","PosCont2") # Deze kolommen willen we niet meenemen
)
colnames(od.long) <- c("Stam", "Herhaling", "conc", "OD") # Nodig om variabele "times"te hernoemen naar "conc"
head(od.long)## Stam Herhaling conc OD
## 1.1.10 1 1 10 0.5676
## 1.2.10 1 2 10 0.5721
## 1.3.10 1 3 10 0.0731
## 1.4.10 1 4 10 0.7895
## 2.1.10 2 1 10 0.8288
## 2.2.10 2 2 10 0.4495Nu kunnen we een plot maken per stam van OD vs. concentratie. We doen hier stam 1.
st1 <- od.long[od.long$Stam==1,]
st1 #check## Stam Herhaling conc OD
## 1.1.10 1 1 10 0.5676
## 1.2.10 1 2 10 0.5721
## 1.3.10 1 3 10 0.0731
## 1.4.10 1 4 10 0.7895
## 1.1.50 1 1 50 0.2304
## 1.2.50 1 2 50 0.3814
## 1.3.50 1 3 50 0.0549
## 1.4.50 1 4 50 0.2330
## 1.1.100 1 1 100 0.1103
## 1.2.100 1 2 100 0.3052
## 1.3.100 1 3 100 0.0526
## 1.4.100 1 4 100 0.1397
## 1.1.250 1 1 250 0.0566
## 1.2.250 1 2 250 0.0673
## 1.3.250 1 3 250 0.0569
## 1.4.250 1 4 250 0.0613
## 1.1.500 1 1 500 0.0499
## 1.2.500 1 2 500 0.0520
## 1.3.500 1 3 500 0.0997
## 1.4.500 1 4 500 0.0583plot(
st1$OD ~ st1$conc, # y ~ x
pch = 19,
main = "Stam 1",
ylab = "OD",
xlab = "concentratie AB"
)
Vraag 14. Beschrijf de relatie tussen OD en concentratie van het antibioticum. Lijkt het alsof stam 1 gevoelig is voor het AB?
Als laatste combineren we alle gegevens. We voegen lijnen toe aan het plaatje voor het gemiddelde en 95% BHI van de positieve en van de negatieve controles. Waarom moeten we nu de grenzen van de y-as veranderen?
# Basis plot met punten
plot(
st1$OD ~ st1$conc, # y ~ x
ylim = c(Negll.nc - 0.1, PosS1ul.nc + 0.1),
pch = 19,
main = "Stam 1",
ylab = "OD",
xlab = "concentratie AB"
)
# Lijnen voor positieve controles
abline(h=PosS1mean.nc, col="green") # Mean
abline(h=PosS1ll.nc, col="green", lty=2) # Lower limit
abline(h=PosS1ul.nc, col="green", lty=2) # Upper limit
# Lijnen voor negatieve controles
abline(h=Negmean.nc, col="red")
abline(h=Negll.nc, col="red", lty=2)
abline(h=Negul.nc, col="red", lty=2)
Vraag 15. Als we groeiremming definiëren als ‘OD lager dan de bovengrens van het 95% BHI van de negatieve controles’, wat moeten we concluderen bij stam 1?
Vraag 16. Herhaal het plaatje voor stammen 2 en 3 en trek conclusies.
Vraag 17a. Stel dat we de negatieve controles zouden gebruiken van alle studenten die het practicum doen. Hoe zal dit de grenzen van het 95% BHI beïnvloeden?
Vraag 17b. Er zijn twee belangrijke aannames voor het berekenen van een BHI. Hoe denk je dat het zit met deze twee aannames met de negatieve controles wanneer we data gebruiken van alle studenten?