R (język programowania)

Ta strona jest wstępem do programowania w R , czyli języku programowania i środowisku do obliczeń statystycznych i wizualizacji wyników.

R-Logo

Poniższe przykłady pokazują pracę w konsoli.

Instalacja

Ubuntu

sudo apt install r-base r-base-dev -y

lub

sudo apt install r-cran-littler

Wprowadzenie

Przykładowa sesja

Rozpoczęcie pracy z konsoli :

R

zakończenie pracy :

q()

co jest skrótem od :

quit()

Możemy od razu odpowiedzieć na pytanie czy chcemy zachować sesję

quit(save="n")

lub

q(save="n")

Edycja

domyślny edytor tekstowy w konsoli

Sprawdzamy jakie jest domyślne ustawienie :

getOption('editor')

przykładowy wynik :

[1] "vi"

czyli w tym przypadki Vim, aby wyjść z Vima wprowadź

:q

Uruchomiamy domyślny edytor tekstowy :

edit()

lub

fix()

albo jednoznacznie :

vi(name = NULL, file = "") 

emacs(name = NULL, file = "") 

pico(name = NULL, file = "") 

xemacs(name = NULL, file = "") 

xedit(name = NULL, file = "")

Inne edytory

• RStudio
• Tinn-R
• RGedit

• Emacs with ESS (“Emacs Speaks Statistics”)
• RGui

Dane

typy danych

dane ucięte

Funkcje

Funkcje :

• data
• head
• names
• dim
• attach i detach

Źródła danych

Dane wewnętrzny R

Pakiet R zawiera przykładowe zbiory danych.

Wyświetlamy listę dostępnych zbiorów danych :

data()

oraz listę zbiorów zawartych we wszystkich dostępnych bibliotekach :

data(package = .packages(all.available = TRUE))

aml

Dane są zawarte w pakiecie survival. Z uwagi na użycie lazyData

nie wczytujemy ich za pomocą data tylko :

aml

cats

Wybieramy zbiór cats z pakietu MASS ( ładujemy pakiet i wczytujemy zbiór ) :

library(MASS)
data(cats)

Sprawdzamy nagłówki :

> names(cats)
[1] "Sex" "Bwt" "Hwt"

Podsumowanie :

 > summary(cats)
 Sex         Bwt             Hwt       
 F:47   Min.   :2.000   Min.   : 6.30  
 M:97   1st Qu.:2.300   1st Qu.: 8.95  
        Median :2.700   Median :10.10  
        Mean   :2.724   Mean   :10.63  
        3rd Qu.:3.025   3rd Qu.:12.12  
        Max.   :3.900   Max.   :20.50

mtcars

Często używany zbiór danych mtcars znajduje się w bibliotece datasets. Nie zawsze jest ona domyślnie ładowana. Wtedy ładujemy bibliotekę :

library(datasets)

i wczytujemy zbiór mtcars :

data(mtcars)

i już możemy go używać. Jeśli chcemy używać nagłówków bez nazwy zbioru danych  :

attach(mtcars)

Funkcji attach() należy jednak używać ostrożnie, ponieważ R w pierwszej kolejności szuka zmiennych w przestrzeni nazw stworzonej poleceniem attach. Gdy chcemy odłączyć dostęp do nazw zmiennych w danym zbiorze stosujemy polecenie detach(). W tym przypadku :

detach(mtcars)

volcano

Informacje topograficzne dotyczące wulkanu Maunga Whau w Auckland: zestaw danych podaje wysokość wulkanu w siatce o wymiarach 10 x 10 metrów. Ten zbiór to macierz z 61 kolumnami i 87 wierszami, w której przechowywane są liczby podwójnej precyzji ( ang. doubles)

Dane w plikach

Dane z plików wczytujemy i zapisujemy do zmiennych. Typ zmiennej zależy od struktury danych. Najłatwiejsze jest wczytywanie danych z plików tekstowych, jak csv. Jest to stosowany przy małej i średniej ilości danych. Przy dużej ilości danych szybsze i objętościowo mniejsze będą pliki binarne.

Katalogi

Najpierw sprawdźmy jaki jest aktualny katalog roboczy :

getwd()

Katalog roboczy możemy zmienić na inny za pomocą :

setwd("nazwa_katalogu")

Potem sprawdźmy jakie są tam pliki :

list.files()

Pliki

CSV

Pliku csv jest plikiem tekstowym.

Tworzenie :

Ręczne

Tworzymy wektor, wyświetlamy go i zapisujemy do pliku :

> wektor=c(1:30)
> wektor
[1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[26] 26 27 28 29 30
> cat(wektor, file = "w.csv",append=F, sep=",")

Tworzymy macierz, wyświetlamy ją i zapisujemy do pliku :

> macierz=matrix(1:6,2,3)
> macierz
[,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> write.table(macierz,file="m.csv", sep=",")

Tworzymy ramkę danych i zapisujemy ją do pliku

> v1 = 1:5
> v2 = c(T,T,F,F,T)
> df = data.frame(v1,v2)
> write.csv(df, file="df.csv")

> df = data.frame(foo=1:5,bar=c(T,T,F,F,T))
> write.csv(df, file="df.csv")

MySQL MySQL pozwala eksportować tabele do pliku csv

W edytorze tekstów

Plik csv możemy utworzyć w każdym edytorze tekstów. Zapisujemy go jak o plik tekstowy z rozszerzeniem csv i przecinkiem jako znakiem rozdzielającym pola.

Konwersja Plik :

• baz danych, np. Accessa z rozszerzeniem mdb
• arkuszy kalkulacyjnych, np. Excela

możemy przekształcić na format csv.

W Accesie i Excelu po utworzeniu bazy danych lub arkusza wybieramy opcję eksportuj do pliku tekstowego. Możemy :

• zmienić znak oddzielający pola ( w polskiej wersji jest średnika a angielskiej przecinek )
• wstawić do pliku nagłówki.

Plik csv z przecinkami nie wczytamy ponownie do Accessa ( błąd), dlatego bezpieczniej jest używać pliku ze średnikiem jako znakiem oddzielającym pola

Otwieranie

Plik lokalny Sposób otwarcia pliku zależy od jego zawartości.

Pliku csv może zawierać dane :

• tego samego typu ułożone w formie tabeli, co w R nazywane jest macierzą ( ang. table)
• różnego typu ułożone w formie tabeli, co w R nazywane jest ramką danych ( ang. data frame)

W pakiecie base jest funkcja do otwarcia pliku csv zawierającego macierz :

read.table("nazwa_pliku") 

W pakiecie utils jest funkcja do otwarcia pliku csv zawierającego ramkę danych z przecinkiem:

read.csv("nazwa_pliku") 

i średnikiem jako znakiem oddzielającym pola :

read.csv2("nazwa_pliku") 

Przykłady użycia :

macierz = read.table("m.csv")
ramka = read.csv("r.csv")

Plik na innym komputerze

Przemysław Biecek udostępnia pliki przykładowymi danymi. Wybieramy plik daneSoc.csv, wczytujemy i zapisujemy do zmiennej ramka :

ramka = read.csv("http://www.biecek.pl/R/dane/daneSoc.csv", sep=";")

Wyświetlamy zawartość zmiennej :

ramka

Inne źródła danych :

• UC Irvine Machine Learning Repository
• datasets at Department of Informatics Nicolaus Copernicus University

Dane tworzone przez użytkownika

> macierz=matrix(1:6,2,3)
> macierz
[,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

Eksport danych

Zapis danych do pliku csv

t=table(a$ROK_WYPISU,a$ICD9CM);
t;
write.table(t,"t.csv", quote = FALSE, sep = ";");

Przygotowanie danych do analizy

Obiekty

Sprawdzamy objekty :

objects()

przykłądowy wynik :

character(0)

Typ obiektu

Sprawdzamy jakiego typu jest zbiór :

> class(cats)
[1] "data.frame"

Jest to ważne ponieważ funkcje odpowiednie dla ramki danych nie muszą być odpowiednie dla macierzy i odwrotnie.

Istnieją funkcje, które można wykonać na obiektach dowolnego typu :

• fix(obiekt)
• class(obiekt)
• unclass(obiekt)
• typeof(obiekt)
• mode(obiekt)
• attributes(obiekt)

Przydatne funkcje :

• names() podaje nazwę zmiennej
• class() podaje typ zmiennej
• str() pokazuje strukturę zbioru danych
• summary() podaje średnią,
• describe() z pakietu Hmisc podaje więcej szczegółów niż summary()
• contents() z pakietu Hmisc
• dims() z pakietu Zelig
• descr() z pakietu descr daje :
  • min, max, mean i quartiles dla zmiennych ciągłych,
  • tablice częstości dla zmiennych czynnikowych
  • długość dla wektorów znaków
• whatis() z pakietu YaleToolkit daje opis zbioru danych

Przykłady pracy z danymi

Wybieramy zbiór cats z pakietu MASS.

Znajdujemy podzbiór kotów płci żeńskiej :

b = subset(cats,Sex=='F')

Liczba kotów płci żeńskiej :

nrow(b)

Otrzymujemy :

47

Tworzymy nową tabelę pokazującą zależność płci od wagi ciała ( Bwt) :

d = table(cats$Sex, cats$Bwt)

Otrzymujemy :

    2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9  3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9
 F  3   9   6  12   4   2   3   3   0   3  2   0   0   0   0   0   0   0   0   0
 M  2   1   8   1   5   8   6   9   7   5  9   6   6   5   5   5   4   1   2   2

Znajdujemy liczbę kotów płci męskiej o największej masie ciała :

d[2,ncol(d)]

Otrzymujemy :

2

Koty płci żeńskiej są w pierwszym rzędzie. Sprawdzamy ich liczbę :

 sum(d[1,])

Otrzymujemy

47

Obliczamy liczbę kotów o najmniejszej masie ciała ( są w 1-szej kolumnie) :

sum(d[,1])

Analiza statystyczna danych

Statystyki opisowe

Średnia masa ciała wszystkich kotów w kg :

> mean(cats$Bwt)
[1] 2.723611

Średnia masa serca wszystkich kotów w g :

> mean(cats$Hwt)
[1] 10.63056

To samo w inny sposób :

> summary(cats)
Sex         Bwt             Hwt       
F:47   Min.   :2.000   Min.   : 6.30  
M:97   1st Qu.:2.300   1st Qu.: 8.95  
       Median :2.700   Median :10.10  
       Mean   :2.724   Mean   :10.63  
       3rd Qu.:3.025   3rd Qu.:12.12  
       Max.   :3.900   Max.   :20.50

Średnia masa ciała wszystkich kotów płci żeńskiej w kg :

mean(subset(cats$Bwt,Sex=='F'))
[1] 2.359574

Analiza przeżyć

Grafika - prezentacja wyników analizy danych

Wczytujemy przykładowe dane :

r = read.csv("http://www.biecek.pl/R/dane/daneSoc.csv", sep=";")

Sprawdzamy nagłówki ramki :

names(r)
[1] "wiek"                  "wyksztalcenie"         "st.cywilny"           
[4] "plec"                  "praca"                 "cisnienie.skurczowe"  
[7] "cisnienie.rozkurczowe"

Możemy tworzyć różne wykresy :

Histogram

Rysujemy histogram wieku :

 hist(r$wiek)

i histogram wieku z krzywą rozkładu normalnego :

 x <- r$wiek
h<-hist(x, breaks=10, col="red", xlab="Wiek",
  main="Histogram with Normal Curve")
xfit<-seq(min(x),max(x),length=40)
yfit<-dnorm(xfit,mean=mean(x),sd=sd(x))
yfit <- yfit*diff(h$mids[1:2])*length(x)
lines(xfit, yfit, col="blue", lwd=2)

Wykres pudełkowy

Rysujemy zależność wieku od stanu cywilnego ( zmiennej liczbowej od czynnikowej) :

boxplot(r$wiek ~ r$st.cywilny)

lub wieku od wykształcenia (dla wszystkich typów) :

plot(wiek~wyksztalcenie, data=r)

lub wieku od wykształcenia podstawowego :

plot(wiek~wyksztalcenie, data=r, subset=wyksztalcenie=="podstawowe")

Wykres workowy

Jest to rozszserzenie do 2 wymiarów wykresu pudełkowego :

library(aplpack)
attach(r)
bagplot(cisnienie.rozkurczowe,cisnienie.skurczowe, xlab="cisnienie.rozkurczowe", ylab="cisnienie.skurczowe",  main="Przyklad bagplot ")

Wykres skrzypcowy

install.packages("UsingR")
library(UsingR)
simple.violinplot(wiek ~wyksztalcenie, data=r, col="red")

Wykres warunkowy

( ang. Conditioning Plot = coplot):

coplot(wiek ~ cisnienie.skurczowe|plec, data = r,   panel = panel.smooth, rows = 1)

Wykres par

pairs(r)

Wykres gęstości

install.packages("UsingR")
library(UsingR)
simple.densityplot(cisnienie.rozkurczowe~wyksztalcenie, data = r)

Wykres

Sprawdzamy jak wygląda zależność pomiędzy ciśnieniem skurczowym i wiekiem ( dwoma zmiennymi liczbowymi) :

plot(r$wiek ~ r$cisnienie.skurczowe)

gradient koloru

Gradient koloru ( ang. colomap)

pals

Pakiet pals

Zależności:

• rgl

Mapy kolorów w pakiecie pals są zdefiniowane:

• za pomocą wzorów matematycznych (np. mapa kolorów Cubehelix), tylko kilka
• większość map kolorów jest pierwotnie zdefiniowana jako gładka krzywa przechodząca przez sekwencję 256 kolorów.

# Install the released version from CRAN:
install.packages("pals")
# Loading required package: pals
require(pals) 
# Show palettes and colormaps as colored bands of below colormaps
pal.bands(coolwarm, parula, ocean.haline, brewer.blues, cubicl, kovesi.rainbow, ocean.phase, brewer.paired(12), stepped, brewer.seqseq2,main="Colormap suggestions")
# increase number of bands to 200
pal.bands(coolwarm, viridis, parula, n=200)
# Show the amount of red, green, blue, and gray in colors of a palette. The gray line corresponds to luminosity.
pal.channels(parula, main="parula")
# Pokaż mapę kolorów za pomocą wykresu czułości kontrastu Campbella-Robsona ( ang. Campbell-Robson Contrast Sensitivity Chart
pal.csf(parula, main="parula")
# The palette is converted to RGB or LUV coordinates and plotted in a three-dimensional scatterplot. The LUV space is probably better, but it is easier to tweak colors by hand in RGB space.
pal.cube(cubehelix)
pal.cube(coolwarm)
# Display a palette on a choropleth map similar to the ColorBrewer website.
pal.map(brewer.paired, n=12, main="brewer.paired")
# A single palette/colormap is shown as five colored bands: 1.Without any modifications. 2. As seen in black-and-white as if photocopied. 3. As seen by the deutan color-blind. 4. As seen by the protan color-blind. 5. As seen by the tritan color-blind.
pal.safe(parula, main="parula")
# compare 3 colormaps
op <- par(mfrow=c(3,1), mar=c(1,1,2,1))
pal.sineramp(jet, main="jet")
pal.sineramp(tol.rainbow, main="tol.rainbow")
pal.sineramp(kovesi.rainbow, main="kovesi.rainbow")
# test colormap
pal.test(parula)
pal.test(viridis) # dark colors are poor
# rysuje zbiór danych volcano z użyciem mapy kolorów parula
pal.volcano(parula)
# Catalog of colormaps and palettes

Pakiety

• r-colorspace i opis
• cetcolor = CET Perceptually Uniform Colour Maps (cetcolor)

rgl

Pakiet rgl to system wizualizacji 3D oparty na OpenGL. Zapewnia interfejs średniego i wysokiego poziomu do użytku w języku R, obecnie wzorowany na klasycznej grafice R, z rozszerzeniami umożliwiającymi interakcję.

install.packages("rgl")
require(rgl)

Pomoc

Gdzie szukać pomocy

Dokumentacja

Sprawdzamy informacje o funkcji "is" :

??is

Otrzymujemy okno z wynikami szukania. Z okna wychodzimy Shift+q

Możemy też użyć funkcji help :

help(is)

kończymy za pomocą klawisza q

Aby się znaleźć w pomocy informacje związane z csv wpisujemy w konsoli :

??csv

Otrzymujemy :

Help files with alias or concept or title matching 'csv' using regular
expression matching:
Hmisc::dataDensityString
                        Internal Hmisc functions
Hmisc::csv.get          Read Comma-Separated Text Data Files
MCMCpack::MCMCSVDreg    Markov Chain Monte Carlo for SVD Regression
survey::dimnames.DBIsvydesign
                        Dimensions of survey designs
survey::open.DBIsvydesign
                        Open and close DBI connections
survey::update.survey.design
                        Add variables to a survey design
utils::read.table       Data Input
utils::write.table      Data Output


Type '?PKG::FOO' to inspect entry 'PKG::FOO TITLE'.

Konwencje nazw

W dokumentacji :

vioplot {vioplot}

oznacza funkcję vioplot z pakietu/biblioteki vioplot.

Biblioteki i pakiety

Przy starcie programu jest wyświetlana informacja jakie biblioteki zostały załadowane. Potem możemy to sprawdzić za pomocą :

sessionInfo()

Sprawdzamy zainstalowane biblioteki :

library()

Wynik jest wyświetlany w domyślnym edytorze

Zainsatlowane pakiety :

installed.packages()

W dokumentacji :

vioplot {vioplot}

oznacza funkcję vioplot z pakietu/biblioteki vioplot. Jeśli ta biblioteka nie jest instalowana przy starcie to należy :

• zainstalować pakiet vioplot :

install.packages("vioplot")

• załadować do pamięci bibliotekę vioplot :

library(vioplot)

• już można używać funkcji vioplot :

vioplot(r$wiek[r$wyksztalcenie=="podstawowe"], r$wiek[r$wyksztalcenie=="srednie"])

Lub za pomocą :

if(!require(somepackage)){
   install.packages("somepackage")
   library(somepackage)

}

obiekty

Funkcja

ls()

lub

objects()

wyświetla zdefiniowane obiekty

Internet

• Forum Miłośników Statystyki ( w tym i R )
• Wrocławski Zlot Użytkowników R
• wyszukiwarka R
• R code at google

Jak zadawać pytania

Jest to opisane w kilku poradnikach.

Dane

W pytaniu podaj dane których używasz. Dzięki temu inni użytkownicy będą mogli sami przetworzyć te dane i porównać wyniki.

• Mogą być to dane tworzone losowo, np. z użyciem funkcji set.seed. Oczywiście musisz podać kod za pomocą którego utworzono dane.
• duży zbiór danych możesz umieścić w przestrzeni publicznej, np. Dropbox
• możesz podać tekstowy dokładny opis danych na których pracujesz utworzony za pomocą funkcji dput. Opis ten można użyć do odtworzenia danych ( dget z pliku albo wkleić tekst structure(....) i przypisać go do zmiennej.

Przykład tekstowego opisu danych:

a =  structure(list(id = c(39L, 471L, 497L, 927L, 663L, 525L, 580L, 
622L, 48L, 727L), age = c(22.6273628946641, 35.237619316895, 
29.6406734238401, 49.7885287820185, 42.6482063541433, 35.8383991257624, 
33.1517001030015, 36.814031442543, 42.5628727298572, 48.4957262906764
)), .Names = c("id", "age"), row.names = c(39L, 471L, 497L, 927L, 
663L, 525L, 580L, 622L, 48L, 727L), class = "data.frame")

Kod

Pokaż kod na którym pracujesz. Podaj minimalną jego ilość, wystarczającą do pokazania błędu/problemu którego nie potrafisz rozwiązać. Kod powinien pracować po skopiowaniu i wklejeniu na czystej sesji R. Nie zapomnij podać niezbędnych bibliotek.

Słownictwo

Używaj właściwego słownictwa. Np jeśli nie wiesz jakiego typu jest obiekt to sprawdź to za pomocą funkcji class.

Dodatkowe materiały

• Portal Edukacyjny Statystyki
• Statystyka matematyczna
• Fizyka statystyczna