Pakiet R w analizie uk³adów z³o¿onych

LABORATORIUM 4

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 kalwiszem ENTER. Wprowadzanie koñczy siê wci¶niêciem klaswisza ENTER w nowej linii. Mo¿na równie¿ przekierowaæ tak wpisywany strumieñ danych do zmiennej.

# PRZYK£AD 4.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 4.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ê¶ciej 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 "test1.dat"

1 "Patient 1" 2 0.5
20 "Patient 10" 10 0.11111
30 "No name" 1 0.99
# PRZYK£AD 4.3

df <- read.table("test1.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 4.4

df <- read.table("test1.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± opcj± 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.

# plik "test2.dat"

name degree cluster
1 "Patient 1" 2 0.5
20 "Patient 10" 10 0.11111
30 "No name" 1 0.99
# PRZYK£AD 4.4

df <- read.table("test2.dat")
df
##          name degree cluster
## 1   Patient 1      2 0.50000
## 20 Patient 10     10 0.11111
## 30    No name      1 0.99000
colnames(df)
## [1] "name"    "degree"  "cluster"
rownames(df)
## [1] "1"  "20" "30"

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

# PRZYK£AD 4.5

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 4.6

df <- read.table("http://www.if.pw.edu.pl/~julas/CSAR/test2.dat", stringsAsFactors = FALSE)
df$name
## [1] "Patient 1"  "Patient 10" "No name"
typeof(df$name)
## [1] "character"

Wczytywanie z bazy danych

Os³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 sie z baz±, a nastêpnie 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 4.7

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 4.8

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 bezposrednia metoda, polegaj±ca na zapisaniu zmiennej do pliku poprzez instrukcjê save(). W ten sposób nie trzeba siê zastanawiaæ nad formatem zapisu.

# PRZYK£AD 4.9

save(dff, file="df")
ls()
## [1] "A"     "con"   "data"  "data1" "df"    "dff"   "pass"  "q"     "user"
rm(dff)
ls()
## [1] "A"     "con"   "data"  "data1" "df"    "pass"  "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 4.10

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.162712  0.256771  3.932454  5.208840  9.115387 11.107331 11.052426
##  [8] 14.531997 15.378468 18.329626 19.151421 22.164435 22.615755 25.958672
## [15] 29.453169 28.215535 32.846416 32.908781 35.030878 39.525842
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ó³cznniki i wyre¶liæ prost± dopasowania.

# PRZYK£AD 4.11

xy.lm <- lm(y ~ x)
summary(xy.lm)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7282 -0.7475 -0.1421  0.8369  1.5390 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.52944    0.49891   -5.07 7.98e-05 ***
## x            2.02957    0.04165   48.73  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.074 on 18 degrees of freedom
## Multiple R-squared:  0.9925, Adjusted R-squared:  0.9921 
## F-statistic:  2375 on 1 and 18 DF,  p-value: < 2.2e-16
xy.lm$coefficients
## (Intercept)           x 
##   -2.529442    2.029573
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 4.12

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 
## -2.4392 -1.4091 -0.5161  1.3540  3.4888 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.4702     1.2208   0.385 0.710183    
## log(x)       -1.8401     0.2976  -6.183 0.000264 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.083 on 8 degrees of freedom
## Multiple R-squared:  0.8269, Adjusted R-squared:  0.8053 
## F-statistic: 38.23 on 1 and 8 DF,  p-value: 0.0002643
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 4.13

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.000000 10.000000 10.000000
## mean           12.866834 16.617981 19.238234 33.420265 44.659542 59.451626
## Std.Dev.        8.052346 11.113669 11.677907  8.765358  9.309627  9.456821
## min            -5.214877  2.050519  7.395381 20.323558 34.413463 39.866017
## Q1             11.307590  9.031966 10.688848 30.060246 36.954333 56.124337
## median         15.628031 15.645531 16.166974 35.080583 42.444097 57.636065
## Q3             18.074671 25.055784 23.030244 37.207205 51.992984 66.635810
## max            21.002211 32.779223 44.146984 50.085231 59.543690 71.270535
## missing values  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000
##                       7        8         9        10
## N              10.00000 10.00000  10.00000  10.00000
## mean           66.88769 79.04111  87.08562  93.70237
## Std.Dev.       15.12287 10.53471  13.83586  11.66258
## min            47.41790 67.60398  63.44521  71.48082
## Q1             54.01120 71.41284  77.95388  92.47812
## median         65.23817 76.34103  89.19307  96.02512
## Q3             79.47495 84.47173  96.97737  99.56108
## max            87.73140 97.64469 107.38068 111.23758
## missing values  0.00000  0.00000   0.00000   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³edów - mo¿na to wykonaæ za pomoc± funkcji errbar() z biblioteki Hmisc.

# PRZYK£AD 4.14

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,]
errbar(x.sb, y.sb, y.sb+e.sb/sqrt(N.sb), y.sb-e.sb/sqrt(N.sb))
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)

Dystrybuanta empiryczna

Jednym z pierwszych wykonywanych wykresów w analizie danych jest porównanie danych z (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 do³adniej przeci±¿onej funkcji plot() na takim obiekcie.

# PRZYK£AD 4.20

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 kodcie wykorzystamy funkcjê par() - s³u¿y ona do przekazywania dodatkowych parametrów do funkcji plot(). W tym przypadku stworzymy ramkê 2 x 2 (mfrow = c(2,2)), która da nam szanse wstawienia kolejnych paneli na rysunek.

# PRZYK£AD 4.21

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)

invisible(sapply(N, make.plots))