Zastosowanie pakietu R w statystyce medycznej

LABORATORIUM 2

Funkcje obslugi macierzy

Część funkcji poznanych wcześniej przy obsłudze wektorów można zastosować także do macierzy, przy czym w przypadku which.min() oraz which.max() konieczne jest użycie dodatkowo funkcji arrayInd() do określenia indeksów macierzy - w przeciwnym wypadku otrzymamy tylko indeks "wektorowy".

# Przykład 3.1
A <- matrix(1:16, 4, 4)
A
##      [,1] [,2] [,3] [,4]
## [1,]    1    5    9   13
## [2,]    2    6   10   14
## [3,]    3    7   11   15
## [4,]    4    8   12   16
sum(A)
## [1] 136
mean(A)
## [1] 8.5
sd(A)
## [1] 4.760952
min(A)
## [1] 1
max(A)
## [1] 16
which.min(A)
## [1] 1
which.max(A)
## [1] 16
arrayInd(which.min(A), dim(A))
##      [,1] [,2]
## [1,]    1    1
arrayInd(which.max(A), dim(A))
##      [,1] [,2]
## [1,]    4    4

Funkcja APPLY

Aby wyznaczyć brzegowe wartości dla macierzy (np. sumę, średnią etc) wykorzystuje się funkcje apply(macierz,liczba,funkcja), przy czym liczba określa, czy odnosimy się do wiersza (1) czy kolumny (2).

# Przykład 3.2
apply(A, 1, sum)
## [1] 28 32 36 40
apply(A, 2, sum)
## [1] 10 26 42 58
apply(A, 1, mean)
## [1]  7  8  9 10
apply(A, 2, sd)
## [1] 1.290994 1.290994 1.290994 1.290994

Zamiast wbudowanych funkcji można także wykorzystac własne funkcje

# Przykład 3.3

# Wlasna funkcja
f <- function(x) {
  sum(x^2)
}

apply(A, 1, f)
## [1] 276 336 404 480
apply(A, 2, f)
## [1]  30 174 446 846

Funkcja SAPPLY

Jak już zostało wcześniej wspomniane, w R należy unikać stosowania pętli. Paradygmat jest następujący: wszędzie, gdzie się da, należy wykorzystywać wektory (lub listy, macierze, ramki danych) i na tych strukturach dokonywać określonych operacji, dostając na wyjściu znów wektor etc. W najprostszym przypadku oznacza to po prostu użycie znanej funkcji (np. sin()) do wszystkich elementów wektora:

# Przykład 3.4

x <- 10
sin(x)
## [1] -0.5440211

Powyższy przykład można również zapisac za pomocą funkcji sapply, która jako argumenty przyjmuje wektor oraz funkcje, która zostanie zastosowana element po elemencie:

# Przykład 3.5

x <- 10
sapply(x, sin)
## [1] -0.5440211

Oczywiście, w tym konkretnym przypadku wygodniej jest użyć opcji z Przykładu 3.4. Gdy jednak sami chcemy skonstruowac funkcje, musimy uciec się do sapply, w przeciwnym wypadku R będzie chciał wykonać naszą funkcję tylko dla pierwszego elementu .

# Przykład 3.6a
# Własna funkcja
g <- function(x) {
  
  y <- sum(1:x)
  if(y > 2 * x) return(x)
  else return(0)
  
}

# Glowny kod
x <- 1:10

# Bezpośrednie użycie g do wektora
g(x)
## Warning in 1:x: numerical expression has 10 elements: only the first used
## Warning in if (y > 2 * x) return(x) else return(0): the condition has
## length > 1 and only the first element will be used
## [1] 0
# Użycie sapply
sapply(x, g)
##  [1]  0  0  0  4  5  6  7  8  9 10

Jeszcze lepiej jest użyć funkcji sapply wewnątrz swojej własnej funkcji

# Przykład 3.6b
# Wlasna funkcja
g <- function(x) {
  
  y <- sapply(x, function(k) sum(1:k))
  ifelse(y > 2 * x, x, 0)
}

# Glowny kod
x <- 1:10

# Bezposrednie użycie g do wektora
g(x)
##  [1]  0  0  0  4  5  6  7  8  9 10

Funkcje anonimowe

Podobnie jak w innych językach skryptowych istnieje możliwość tworzenia funkcji anonimowych. Są to najczęściej krótkie wyrażenia, nie tylko arytmetyczne, bez nazwy (stąd anonimowe), które istnieją ulotnie - po wykonaniu operacji zostają usunięte z pamięci.

# Przykład 3.7
a <- 1:10
sapply(a, function(x) x^2 - 2)
##  [1] -1  2  7 14 23 34 47 62 79 98

Za pomocą funkcji anonimowych oraz zagnieżdżonych funkcji sapply można tworzyc alternatywę dla podwójnych pętli.

# Przykład 3.8

a <- 1:10
sapply(a, function(x) {sapply(rev(a), function(y) x * y)})
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]   10   20   30   40   50   60   70   80   90   100
##  [2,]    9   18   27   36   45   54   63   72   81    90
##  [3,]    8   16   24   32   40   48   56   64   72    80
##  [4,]    7   14   21   28   35   42   49   56   63    70
##  [5,]    6   12   18   24   30   36   42   48   54    60
##  [6,]    5   10   15   20   25   30   35   40   45    50
##  [7,]    4    8   12   16   20   24   28   32   36    40
##  [8,]    3    6    9   12   15   18   21   24   27    30
##  [9,]    2    4    6    8   10   12   14   16   18    20
## [10,]    1    2    3    4    5    6    7    8    9    10

Tworzenie wykresów

Standardowa funkcja wykorzystywana do tworzenia wykresów jest plot(x,y), gdzie x i y sa odpowiednio wektorami liczb.

# Przykład 3.9
x <- 1:10
plot(x, x^2)

Funkcja plot() ma pokaźny zestaw różnych opcji, przy czym najczęściej wykorzystywane to xlab="..." (tytuł osi X), ylab="..." (tytuł osi Y), main="..." (tytuł wykresu), col="..." (kolor punktów), pch=... (kształt punktów), cex="..." (rozmiar punktów), font=... (typ czcionki osi: 1 - normalna, 2 - pogrubiona, 3 - kursywa, 4 - pogrubiona kursywa). W przypadku nazw osi można wykorzystywać funkcje expression(), która koduje wyrażenia matematyczne i litery greckie (np. ^ to indeks gorny, [..] - indeks dolny, itd.).

# Przykład 3.10
x <- 1:10
plot(x, x^2, xlab = "x", ylab = expression(f(xi)==x^2), col = "red", pch = 19, font = 2, font.lab = 4, main = "Wykres funkcji f(x)", font.main = 3, cex = 2)

W przypadku potrzeby zapisania wykresu do pliku należy skorzystać z jednej z komend png(), jpeg() czy tiff(), a następnie zamknąć strumień komendą dev.off().

# Przykład 3.11

png("fig2.png")
plot(x, x^2, xlab = "x", ylab = expression(f(xi)==x^2), col = "red", pch = 19, font = 2, font.lab = 4, main = "Wykres funkcji f(x)", font.main = 3, cex = 2)
dev.off()
## png 
##   2

Do umieszczenia legendy na wykresie służy funkcja legend(), której najważniejszymi argumentami są:

# Przykład 3.12
x <- 1:10
plot(x, x^2, xlab = "x", ylab = expression(f(x)), col = "red", bg="red", pch = 21, cex = 2)
points(x, x^1.5, col="blue", bg="blue", pch=22, cex=2)
legend("topleft", legend=c(expression(x^1.5),expression(x^2)), pch = c(22,19), col = c("blue","red"), pt.bg=c("blue","red"))

Tworzenie histogramów

Najprostszą metodą tworzenia histogramu (w formie tekstowej) jest wykorzystanie funkcji tabulate(). Efektem jej działania jest zliczenie wartości całkowitych zawartych w wektorze. W przypadku liczb rzeczywistych dokonywane jest zaokrąglenie w dół.

# Przykład 3.13

y <- c(0,0,1,2,3,1,2,3,4)
tabulate(y)
## [1] 2 2 2 1
y <- c(0,0,1.1,1.9,2.1,2,3)
tabulate(y)
## [1] 2 2 1

W przypadku histogramu 2D można posłuzyć się funkcją table(), która stworzy dwuwymiarową tablicę współwystępowania wartości.

# Przykład 3.14

df <- data.frame(x=c(1,1,2,2,3,4,5), y=c(2,2,3,1,5,5,5))
df
##   x y
## 1 1 2
## 2 1 2
## 3 2 3
## 4 2 1
## 5 3 5
## 6 4 5
## 7 5 5
table(df)
##    y
## x   1 2 3 5
##   1 0 2 0 0
##   2 1 0 1 0
##   3 0 0 0 1
##   4 0 0 0 1
##   5 0 0 0 1

Jednak klasyczna funkcja odpowiedzialna za tworzenie histogramów jest hist(). Samo jej wywołanie daje efekt narysowania histogramu o zadanej liczbie przedziałów (binów).

# Przykład 3.15
x <- c(1,1,1,2,2,4,10)
hist(x)

W pewnym sensie nawet ważniejszą rzeczą od samego wykresu jest zawartość zmiennej, do której zostanie zapisany jego wynik. Otrzymamy z niej informacje nie tylko o liczbie zliczeń (counts), ale także o granicach binów (breaks), ich środkach (mids) jak również funkcji gęstosci (density).

# Przykład 3.16
x <- c(1,1,1,2,2,4,10)
h <- hist(x)

h
## $breaks
## [1]  0  2  4  6  8 10
## 
## $counts
## [1] 5 1 0 0 1
## 
## $density
## [1] 0.35714286 0.07142857 0.00000000 0.00000000 0.07142857
## 
## $mids
## [1] 1 3 5 7 9
## 
## $xname
## [1] "x"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"

Losowanie wartości

Funkcja SAMPLE

Za pomocą funkcji sample(x) można w pakiecie R uzyskać permutację oryginalnego zbioru (wektora x). W przypadku podania konkretnej liczby próbek, mniejszej niż rozmiar x elementy zostaną wylosowane bez zwracania. Wreszcie, podanie opcji replace=TRUE umożliwi losowanie ze zwracaniem. Dodatkowo, podanie wektora prob daje możliwość sterowania prawdopodobieństwem wylosowania konkretnych elementów wektora x. Funkcja nie ogranicza się jedynie do typów liczbowych.

# Przykład 3.17

sample(1:3, 2)
## [1] 3 1
sample(1:10, 2)
## [1] 5 6
sample(1:10)
##  [1]  8  4  6  2  5  7 10  3  1  9
sample(1:10, 4)
## [1] 10  4  9  2
sample(1:10, 20, replace=TRUE)
##  [1]  8 10  1  3  5 10  9  4  8  6  7  4  8  3  3  1  1  8  2  9
sample(1:3, 10, replace=TRUE, prob=c(0.1,0.8,0.1))
##  [1] 2 2 2 2 3 2 2 2 2 2
sample(letters[1:3], 10, replace=TRUE)
##  [1] "b" "b" "c" "a" "b" "c" "b" "b" "b" "b"
sample(c(0.1, 0.2, 0.3), 10, replace=TRUE)
##  [1] 0.3 0.1 0.2 0.1 0.3 0.3 0.3 0.2 0.2 0.3

Rozkłady prawdopodobieństw

W pakiecie R jest zaimplementowany cały zestaw standardowych rozkładów prawdopodobieństw (rozkład Gaussa, dwumianowy etc), do których odwołujemy się w ten sam sposob r|p|d|qnazwa(), gdzie w miejsce nazwa należy wstawić nazwę odpowiedniego rozkładu prawdopodobieństwa (np. norm, unif, exp, binom), a pierwsza litera oznacza typ funkcji:

  • r - losowanie z rozkładu (np. runif(5) - losowanie pięciu liczb z rozkładu jednostajnego),

  • p - dystrybuantę (np. pnorm(2, 0, 1.5) - wartość dystrybuanty dla rozkładu normalnego o średniej 0 i odchyleniu 1.5 w punkcie 2),

  • d - gęstość prawdopodobieństwa (np. dexp(2, 0.1) - gęstość prawdopodobieństwa dla rozkładu wykładniczego o parameterze 0.1 w punkcie 2),

  • q - kwantyl (np. qnorm(0.95, 0, 1) zwróci kwantyl rzędu 0.95 z rozkładu normalnego o średniej 0 i odchyleniu standardowym 1).

# Przykład 3.18

x <- seq(-2, 2, .1)
plot(x, dnorm(x, 0, 0.5), ylim = c(0,1.5), pch = 19)
points(x, dnorm(x, 0, 1), pch = 19, col = "blue")
points(x, dnorm(x, 0, 0.3), pch = 19, col = "green", t = "o")
lines(x, pnorm(x, 0, 0.3), col = "red", lwd = 2)

# Przykład 3.19

runif(10)
##  [1] 0.6607871 0.7717536 0.7353429 0.5211421 0.3531406 0.9384612 0.9312509
##  [8] 0.6434158 0.3916183 0.7073876
runif(10, 5, 10)
##  [1] 8.269966 6.724866 5.963732 7.468810 8.479985 6.939653 7.321587
##  [8] 8.478122 6.719646 9.895461


Zadanie punktowane

Wylosuj \(n\) liczb z rozkładu normalnego o średniej \(m\) i odchyleniu standardowym \(s\) (\(n\),\(m\),\(s\) powinny być zmiennymi których wartość będzie można zmieniać w skrypcie). Następnie wykonaj histrogram z tych danych i umieść na wykresie empiryczną gęstość prawdopodobieństwa w postaci punktów (tę gęstość można wyciągnąć z obiektu histogram. Dodaj drugą serię do wykresu (czerwona linia), która będzie przedstawiać teoretyczną gęstość prawdopodobieństwa. Podpisz osie i dodaj legendę tak jak na rysunku poniżej.