Zastosowanie pakietu R w statystyce medycznej

LABORATORIUM 3

Instalacja bibliotek

W jêzyku R dodatkowe biblioteki mo¿na zainstalowac bezpo¶rednio w terminalu przy pomocy komendy install.packages("nazwa_biblioteki"), np.

install.packages("fields")
install.packages("Hmisc")

Wbudowane zbiory danych

Po¶ród zbioru podstawowych bibliotek R znajduje siê pakiet datasets, który zawiera kilkadziesi±t gotowych do wykorzystania zbiorów danych. Komenda data() zwraca listê tych zbiorów, natomiast wczytanie do pamiêci konkretnych danych odbywa siê poprzez przekazanie nazwy zbioru jako argument funkcji data("nazwa_zbioru").

# Wypisanie zbiorów z pakietu datasets
data()

# Wczytanie zbioru
data("ChickWeight")

# Wypisanie wszystkich zbiorów ¿e wszystkich dostêpnych pakietów
data(package = .packages(all.available = TRUE))

Wczytywanie danych

Wczytywanie danych z konsoli

Za pomoc± funkcji scan() mo¿na w pakiecie R wczytywaæ dane bezpo¶rednio z konsoli. Kolejne elementy oddzielamy spacjami lub klawiszem Enter. Wprowadzanie koñczymy wci¶niêciem klawisza Enter w nowej linii. mo¿na równie¿ przekierowaæ tak wpisywany strumieñ danych do zmiennej.

# Przyk³ad 3.1

> scan()
1: 5 3
3: 7 8 9 10
7: 
Read 6 items
[1]  5  3  7  8  9 10
> x <- scan()
1: 6 3 9 10 23 -9
7: 
Read 6 items
> x
[1]  6  3  9 10 23 -9

Obs³ugiwane jest tak¿e wprowadzanie innych typów zmiennych poprzez podanie opcji scan(what="...")

# Przyk³ad 3.2

> x <- scan(what="character")
1: a b c g 
5: cos
6: 
Read 5 items
> x
[1] "a"   "b"   "c"   "g"   "cos"

Wczytywanie danych z pliku

Jedn± z najczêsciej wykorzystywanych funkcji do wczytywania danych z pliku jest read.table(), tworz±ca z wczytanego zbioru ramkê danych. Oznacza to, ¿e w pliku ka¿da linia powinna zawieraæ tyle samo pól, a poza tym ka¿da kolumna musi zawieraæ ten sam typ danych.

# plik "ramka1.dat"

1 "Patient 1" 2 0.5
20 "Patient 10" 10 0.11111
30 "No name" 1 0.99
# Przyk³ad 3.3

df <- read.table("data/ramka1.dat")
df
##   V1         V2 V3      V4
## 1  1  Patient 1  2 0.50000
## 2 20 Patient 10 10 0.11111
## 3 30    No name  1 0.99000

W przypadku, gdy chcemy nadaæ nazwy poszczególnym kolumnom ramki, podajemy opcjê read.table(col.names="..."). Spacje w nazwach zostan± zast±pione kropkami.

# Przyk³ad 3.4

df <- read.table("data/ramka1.dat", col.names=c("id", "name", "degree", "clust coeff"))
df
##   id       name degree clust.coeff
## 1  1  Patient 1      2     0.50000
## 2 20 Patient 10     10     0.11111
## 3 30    No name      1     0.99000
df$name
## [1] Patient 1  Patient 10 No name   
## Levels: No name Patient 1 Patient 10

Inn± opcja jest podanie nazw kolumn w samym pliku - je¶li w pierwszym wierszu jest o jedno pole mniej ni¿ w kolejnym, funkcja automatycznie potraktuje pierwsz± liniê jako nazwy kolumn, natomiast pierwsz± kolumnê jako nazwy wierszy. Je¶li w pliku podane s± nazwy kolumn, ale nie ma podanych nazw wierszy (czyli pierwszy wiersz zawiera tê sam± liczbê pól), to nale¿y u¿yæ opcji header=TRUE wewn±trz funkcji read.table(...).

# plik "ramka2.dat"

name degree cluster
1 "Patient 1" 2 0.5
20 "Patient 10" 10 0.11111
30 "No name" 1 0.99

# plik "ramka3.dat"
name age weight
"Kowalski" 38 94.3
"Nowak" 25 67.5
"Malinowski" 49 84.7
# Przyk³ad 3.5

df2 <- read.table("data/ramka2.dat"); df2
##          name degree cluster
## 1   Patient 1      2 0.50000
## 20 Patient 10     10 0.11111
## 30    No name      1 0.99000
colnames(df2)
## [1] "name"    "degree"  "cluster"
rownames(df2)
## [1] "1"  "20" "30"
df3 <- read.table("data/ramka3.dat", header = TRUE); df3
##         name age weight
## 1   Kowalski  38   94.3
## 2      Nowak  25   67.5
## 3 Malinowski  49   84.7

Warto tu zaznaczyæ, ¿e read.table() wczytuje domy¶lnie (niestety) ³añcuchy znaków jako typ wyliczeniowy (factor).

# Przyk³ad 3.6

df$name
## [1] Patient 1  Patient 10 No name   
## Levels: No name Patient 1 Patient 10
typeof(df$name)
## [1] "integer"

Aby temu zaradziæ, nale¿y u¿yæ opcji stringsAsFactors = FALSE. Ponadto, w przypadku gdy znamy adres url pliku, zamiast go pobieraæ i wczytywaæ lokalnie, mo¿na po prostu podaæ jego lokalizacjê:

# Przyk³ad 3.7

df2 <- read.table("http://www.if.pw.edu.pl/~paluch/MSR/data/ramka2.dat", stringsAsFactors = FALSE)
df2$name
## [1] "Patient 1"  "Patient 10" "No name"
typeof(df2$name)
## [1] "character"

Wczytywanie z bazy danych

Obs³ug± bazy danych MySQL zajmuje siê biblioteka RMySQL, któr± wczytujemy komend± library(RMySQL). Dysponuj±c poni¿sz± tabel± w bazie:

# tabela table_test

+----+------------------+--------+
| id | name             | salary |
+----+------------------+--------+
|  1 | Jan Kowalski     |   1500 |
|  2 | Karol Wisniewski |   4000 |
|  3 | Karol Wielki     |   4000 |
+----+------------------+--------+

wykorzystujemy funkcjê dbConnect() do po³±czenia siê z baz±, a nastepniê dbSendQuery() do wys³ania zapytania. Wreszcie za pomoc± komendy fetch() pobieramy rezulaty zapytania do zmiennej typu data frame. Wyczyszczenie rezultatów odbywa siê komend± dbClearResult(), a zamkniêcie po³±czenia z baz± danych dbDisconnect().

# Przyk³ad 3.8

library(RMySQL)
## Loading required package: DBI
user <- "test"; pass <- "CSARuser"
con <- dbConnect(MySQL(), host="194.29.174.45", user=user, password=pass, dbname="CSAR")
q <- dbSendQuery(con, "SELECT * FROM table_test;")
data <- fetch(q)
data
##   id             name salary
## 1  1     Jan Kowalski   1500
## 2  2 Karol Wisniewski   4000
## 3  3     Karol Wielki   4000
is.data.frame(data)
## [1] TRUE
dbClearResult(q)
## [1] TRUE
q <- dbSendQuery(con, "SELECT * FROM table_test WHERE salary > 1500;")
data1 <- fetch(q)
data1
##   id             name salary
## 1  2 Karol Wisniewski   4000
## 2  3     Karol Wielki   4000
dbDisconnect(con)
## Warning: Closing open result sets
## [1] TRUE

Zapisywanie danych

W przypadku macierzy oraz ramek danych odpowiednim sposobem zapisu danych jest u¿ycie funkcji write.table(). W przypadku, gdy dane maj± byæ "czyste" (bez nazw kolumn i wierszy), u¿ywamy opcji write.table(col.names=F, row.names=F).

# Przyk³ad 3.9

dff <- data.frame(x=1:3, names=c("Aaaa", "Bbbb", "Ccc"))
dff
##   x names
## 1 1  Aaaa
## 2 2  Bbbb
## 3 3   Ccc
write.table(dff, "test4.dat")
A <- matrix(1:10, 2, 5)
write.table(A, "test5.dat")
write.table(A, "test6.dat", row.names=F, col.names=F)

Istnieje wreszcie najbardziej bezpo¶rednia metoda, polegaj±ca na zapisaniu zmiennej do pliku poprzez instrukcjê save(). W ten sposób nie trzeba siê zastanawiaæ nad formatem zapisu.

# Przyk³ad 3.10

save(dff, file="df")
ls()
##  [1] "A"     "con"   "data"  "data1" "df"    "df2"   "df3"   "dff"  
##  [9] "pass"  "q"     "user"
rm(dff)
ls()
##  [1] "A"     "con"   "data"  "data1" "df"    "df2"   "df3"   "pass" 
##  [9] "q"     "user"
dff
## Error in eval(expr, envir, enclos): object 'dff' not found
load("df")
dff
##   x names
## 1 1  Aaaa
## 2 2  Bbbb
## 3 3   Ccc

Regresja liniowa

W wielu przypadkach interesuje nas wykonanie analizy regresji opracowywanych danych. W pakiecie R s³u¿y do tego funkcja lm(), przy czym specyficzny jest sposób wprowadzania formu³y - w R wykorzystuje siê symbol tyldy ~ do pokazania zale¿no¶ci pomiêdzy zmiennymi (np. y ~ x oznacza zale¿no¶æ pomiêdzy x i y. Poni¿ej generujemy zaburzon± losowo zale¿no¶æ liniow±:

# Przyk³ad 3.11

x <- 1:20
y <- 2*x-2 + runif(length(x), -3, 3)
x
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
y
##  [1]  1.8192515  0.2375289  5.7683671  6.5007705  5.7581109 11.5472238
##  [7] 11.4347324 14.9502152 15.9868067 18.5926516 19.4592244 22.6583204
## [13] 26.4927403 28.4735556 25.4979993 31.2777868 32.9944871 34.2511302
## [19] 38.3612121 36.7394060
plot(x, y)

Wykonujemy regresjê liniow± i wynik procedury zapisujemy do zmiennej xy.lm. Po wywo³aniu funkcji summary() otrzymamy wszystkie interesuj±ce nas warto¶ci dotycz±ce regresji. Bezpo¶rednie odwo³anie do wspó³czynników otrzymujemy poprzez pole coefficients, przy czym coefficients[1] to punkt przeciêcia, a coefficients[2], to wspó³czynnik kierunkowy. Mo¿emy nastêpnie zapisaæ wspó³czynniki i wykre¶liæ prost± dopasowania.

# Przyk³ad 3.12

xy.lm <- lm(y ~ x)
summary(xy.lm)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0347 -0.9477  0.1909  1.2705  2.0012 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.77612    0.72463  -2.451   0.0247 *  
## x            2.02059    0.06049  33.403   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.56 on 18 degrees of freedom
## Multiple R-squared:  0.9841, Adjusted R-squared:  0.9832 
## F-statistic:  1116 on 1 and 18 DF,  p-value: < 2.2e-16
xy.lm$coefficients
## (Intercept)           x 
##   -1.776115    2.020590
b <- xy.lm$coefficients[1]
a <- xy.lm$coefficients[2]
plot(x, y, pch=19)
lines(x, a*x+b, col="red", lwd=3)

Podobn± procedurê mo¿na przeprowadziæ dla przetransformowanych zmiennych. Tutaj jest to zrandomiznowana zale¿no¶æ \(y \sim x^{-2}\)

# Przyk³ad 3.13

x <- c(1,2,5,10,20,50,100,200,500,1000)
y <- (x*(1+runif(length(x),-0.9,0.9)))^(-2)
plot(x, y, log="xy")

xy.lm <- lm(log(y) ~ log(x))
summary(xy.lm)
## 
## Call:
## lm(formula = log(y) ~ log(x))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4895 -1.1963 -0.6035  0.9406  3.0430 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.2821     0.9799   1.308    0.227    
## log(x)       -2.2317     0.2389  -9.342 1.41e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.672 on 8 degrees of freedom
## Multiple R-squared:  0.916,  Adjusted R-squared:  0.9055 
## F-statistic: 87.28 on 1 and 8 DF,  p-value: 1.407e-05
b <- exp(xy.lm$coefficients[1])
a <- xy.lm$coefficients[2]
plot(x, y, pch=19, log="xy")
lines(x, b*x^a, col="red", lwd=3)

Binowanie danych

W wielu przypadkach konieczne lub po¿±dane jest dokonanie "binowania" danych, czyli podzielenia osi X na przedzia³y, a nastêpnie wyliczenia warto¶ci ¶rednich w tych przedzia³ach. W R mo¿na do tego wykorzystaæ funkcjê stats.bin() z biblioteki fields.

# Przyk³ad 3.14

x <- 1:100
y <- x + runif(length(x), -20, 20)

library(fields)
stats.bin(x, y)
## $centers
##  [1]  5 15 25 35 45 55 65 75 85 95
## 
## $breaks
##  [1]   0  10  20  30  40  50  60  70  80  90 100
## 
## $stats
##                        1         2         3        4        5        6
## N              10.000000 10.000000 10.000000 10.00000 10.00000 10.00000
## mean            6.326604 13.120527 23.586485 29.04415 41.91543 58.92546
## Std.Dev.        8.717828  7.661717 11.593177 12.37924 11.58576 12.28617
## min            -5.813408  4.570535  6.726506 14.70920 24.22519 46.78778
## Q1              1.510255  5.828519 15.902061 20.96305 36.02693 49.77282
## median          3.885198 12.771574 25.228769 27.92984 39.15062 53.98639
## Q3             11.431417 20.626427 29.607007 31.51168 47.77494 69.20096
## max            21.144786 22.214020 44.520810 53.33758 61.72505 79.13718
## missing values  0.000000  0.000000  0.000000  0.00000  0.00000  0.00000
##                       7        8          9        10
## N              10.00000 10.00000  10.000000  10.00000
## mean           65.48497 71.79921  84.331933  94.11269
## Std.Dev.       10.73156 12.46582   8.741132  13.69385
## min            51.16072 53.34288  73.477091  72.46845
## Q1             56.07058 61.54615  77.616545  86.56396
## median         66.84778 70.70682  84.911308  96.72417
## Q3             69.91567 83.98114  88.397075 105.13925
## max            81.92220 85.94273 100.363484 112.04342
## missing values  0.00000  0.00000   0.000000   0.00000
xy.sb <- stats.bin(x, y)
plot(x, y)
with(xy.sb, points(centers, stats[2,], col="red", pch=19, cex=1.5))

Aby mo¿na by³o wykonaæ "przyzwoity" wykres z dopasowaniem liniowym, nale¿a³oby umie¶ciæ na nim tak¿e s³upki b³êdów - mo¿na to wykonaæ za pomoc± funkcji errbar() z biblioteki Hmisc.

# Przyk³ad 3.15

library(Hmisc)
x.sb <- xy.sb$centers
y.sb <- xy.sb$stats[2,]
e.sb <- xy.sb$stats[3,]
N.sb <- xy.sb$stats[1,]
xy.lm <- lm(y.sb ~ x.sb)
a.xy <- xy.lm$coefficients[2]
b.xy <- xy.lm$coefficients[1]
errbar(x.sb, y.sb, y.sb+e.sb/sqrt(N.sb), y.sb-e.sb/sqrt(N.sb))
points(x, y, col="gray", pch=19)
lines(x, a.xy*x+b.xy, col="red", lwd=3, lty=3)

W powy¿szym przyk³adzie funkcja errbar() jest wywo³ana zamiast funkcji plot(). Mo¿liwe jest jednak dodanie s³upków b³êdów jako kolejnej warstwy do utworzonego wykresu na takiej samej zasadzie jak u¿ywamy funkcji points() i lines(). Trzeba pamiêtaæ wtedy o parametrze add = TRUE wewn±trz funkcji errbar(). Kolor s³upków b³êdów mo¿na zmieniæ parametrem errbar.col.

# Przyk³ad 3.16

plot(x, y, col="gray", pch=19, xlab = "x", ylab = "f(x)", main = "Binowanie i regresja liniowa")
errbar(x.sb, y.sb, y.sb+e.sb/sqrt(N.sb), y.sb-e.sb/sqrt(N.sb), col="blue", errbar.col="blue", lwd=1.25, add=TRUE)
lines(x, a.xy*x+b.xy, col="red", lwd=3, lty=3)

Dystrybuanta empiryczna

Jednym z pierwszych wykonywanych wykresów w analizie danych jest porównanie danych (a dok³adniej ich rozk³adu) z jakim¶ konkretnym rozk³adem teoretycznym. Najczê¶ciej jest to dystrybuanta - szybkim sposobem na otrzymanie tego skumulowanego rozk³adu dla danych (czyli dystrybuanty empirycznej) jest wywo³anie funkcji ecdf(), a dok³adniej przeci±¿onej funkcji plot() na takim obiekcie.

# Przyk³ad 3.17

x <- c(1,1,1,2,5,6,3,7,8,10)
plot(ecdf(x))

Przy okazji warto zobaczyæ, jak taka dystrybunata zachowuje siê w przypadku wylosowanych danych o ró¿nym rozmiarze próbki. W poni¿szym kodzie wykorzystamy funkcjê par() - s³u¿y ona do przekazywania dodatkowych parametrów do funkcji plot(). W tym przypadku stworzymy ramkê 2x2 (mfrow = c(2,2)), która da nam szansê wstawienia kolejnych paneli na rysunek.

# Przyk³ad 3.18

get.x <- function(x) {
  return(seq(min(x), max(x), length.out = 100))
}

make.plots <- function(N) {
  x <- rnorm(N, 0, 1)
  plot(ecdf(x), main = N)
  xx <- get.x(x)
  lines(xx, pnorm(xx, 0, 1), col = "red", lwd = 2)  
}

par(mfrow = c(2,2))

N <- c(10, 50, 100, 500)

sapply(N, make.plots)

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL


Zadanie punktowane

Wczytaj do pamiêci zbiór danych ChickWeight i zapoznaj siê z nim, a nastêpnie stwórz ramkê danych zawieraj±c± wy³acznie te rekordy, dla których warto¶æ zmiennej Diet wynosi 1. Wykonaj wykres wagi kur od d³ugo¶ci ich ¿ycia. Zbinuj dane i umie¶æ na wykresie w postaci czerwonych punktów ze s³upkami b³êdów odpowiadaj±cymi niepewno¶ci standardowej (standard error). Dla zbinowanych danych wykonaj regresjê liniow± i wykre¶l prost± najlepszego dopasowania w kolorze niebieskim. Zapisz do pliku tekstowego ramkê danych zawieraj±c± trzy kolumny: wspó³rzêdne x i y punktów oraz niepewno¶æ standardow±. 

##     x        y   std.err
##  1  0 41.40000 0.2224268
##  2  2 47.25000 0.9566251
##  3  4 56.47368 0.9470435
##  4  6 66.78947 1.7796428
##  5  8 79.68421 3.1604765
##  6 10 93.05263 5.1716013