####################################################################################### ################ Statistická laboratoř - 5. workshop ################################## ############### Adéla Vrtková, Martina Litschmannová ################################## ####################################################################################### ## Máme data, a co dál? - 1. část # 1. Spustíme potřebné balíčky, které obsahují další statistické funkce # 2. Nastavíme pracovní adresář, odkud importujeme data, popř. kam chceme ukládat výstupy # 3. Importujeme data (z pracovního adresáře, z internetu) # 4. Pre-processing -> a) Podíváme se na data # b) uložíme si data ve více formátech (každá funkce má "radši" jiný formát) # 5. Analýza kvalitativních proměnných ## Máme data, a co dál? - 2. část # 6. Analýza kvantitativních proměnných # 7. Identifikace a rozhodnutí o vyloučení/ponechání odlehlých pozorování ####################################################################################### ## 3. Načtení datového souboru + 4. Preprocessing #################################### # Načtení a uložení datového souboru ve formátu csv2 z internetu do datového rámce data data=read.csv2(file="http://am-nas.vsb.cz/lit40/DATA/aku.csv") # Uložení prvních 4 sloupců dat. rámce data do dat. rámce data5 data5=data[,c(1:4)] ## Převod dat do standardního datového formátu data5S=reshape(data[,1:4], # část datového rámce, která bude převáděna do std. datového formátu direction="long", # parametr určující tzv. "long" nebo "wide" formát varying=c("A5","B5","C5","D5"), # názvy proměnných, které mají být zařazeny do sloupce hodnot v.names="kap5", # pojmenování sloupce hodnot times=c("A","B","C","D"), # varianty proměnné, která přiřazuje identifikátory jednotlivým hodnotám proměnné kap5 timevar="vyrobce") # pojmenování proměnné obsahující identifikátory proměnné kap5 # převedení proměnné data5S$vyrobce na typ factor data5S$vyrobce=as.factor(data5S$vyrobce) # odstranění nadbytečné proměnné id z datového rámce data5S data5S=data5S[,-3] # odstranění NA z datového rámce data5S data5S=na.omit(data5S) ## Převod párových dat do standardního datového formátu dataS=reshape(data, direction="long", varying=list(c("A5","B5","C5","D5"), c("A100","B100","C100","D100")), v.names=c("kap5","kap100"), times=c("A","B","C","D"), timevar="vyrobce") # převedení proměnné dataS$vyrobce na typ factor dataS$vyrobce=as.factor(dataS$vyrobce) # odstranění nadbytečné proměnné id z datového rámce dataS dataS=dataS[,-4] # odstranění NA z datového rámce dataS dataS=na.omit(dataS) # Definování nové proměnné v stávajícím datovém rámci dataS$pokles=dataS$kap5-dataS$kap100 ## Vytvoření samostatných proměnných a5=dataS$kap5[dataS$vyrobce=="A"] b5=dataS$kap5[dataS$vyrobce=="B"] c5=dataS$kap5[dataS$vyrobce=="C"] d5=dataS$kap5[dataS$vyrobce=="D"] a100=dataS$kap100[dataS$vyrobce=="A"] b100=dataS$kap100[dataS$vyrobce=="B"] c100=dataS$kap100[dataS$vyrobce=="C"] d100=dataS$kap100[dataS$vyrobce=="D"] pokles.a=dataS$pokles[dataS$vyrobce=="A"] pokles.b=dataS$pokles[dataS$vyrobce=="B"] pokles.c=dataS$pokles[dataS$vyrobce=="C"] pokles.d=dataS$pokles[dataS$vyrobce=="D"] ### Poznámky pro zopakování principu grafiky v R ###################################### # základem jsou tzv. high-level funkce, které vytvoří graf (tj. otevřou grafické oknou a vykreslí dle zadaných parametrů) # na ně navazují tzv. low-level funkce, které něco do aktviního grafického okna přidají, samy o sobě neotevřou nové # př. low-level funkcí - např. abline, points, lines, legend, title, axis ... které přidají přímku, body, legendu... # tzn. před použitím "low-level" funkce je potřeba, volat "high-level" funkci (např. plot, boxplot, hist, barplot, pie,...) # další grafické parametry naleznete v nápovědě # nebo např. zde http://www.statmethods.net/advgraphs/parameters.html # nebo zde https://flowingdata.com/2015/03/17/r-cheat-sheet-for-graphical-parameters/ # nebo http://bcb.dfci.harvard.edu/~aedin/courses/BiocDec2011/2.Plotting.pdf ## Barvy v R # http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf # https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/colorPaletteCheatsheet.pdf ####################################################################################### ## 6. Explorační analýza a vizualizace kvantitativní proměnné ######################### ## Popisná statistika ################################################################# summary(dataS$kap5) # Výpočet průměru jedné proměnné mean(dataS$kap5) mean(a5) # Pozor na chybějící hodnoty mean(data$C5) mean(data$C5,na.rm=TRUE) mean(na.omit(data$C5)) # Výpočet mediánu jedné proměnné quantile(dataS$kap5,probs=0.5) quantile(a5,probs=0.5) # Určení rozsahu length(na.omit(dataS$kap5)) ####################################################################################### ## funkce, které umožní aplikaci vybrané funkce na sloupce datového rámce, na hodnoty vektoru apod. x=1:4 # - lapply - na vstup aplikuje zvolenou funkci, vrátí list (seznam) hodnot lapply(x,sqrt) # - sapply - na vstup aplikuje volenou funkci, vrátí vektor hodnot sapply(x,sqrt) # - vapply - stejné jako sapply, akorát s parametrem navíc, který specifikuje výstup vapply(x,sqrt,numeric(1)) # - tapply - lze aplikovat na vektor, jemuž je přiřazen faktor rozlišující hodnoty do skupin # - použijeme dále ####################################################################################### # Výpočet průměrů kapacit po 5 cyklech dle výrobců tapply(data5S$kap5, data5S$vyrobce, mean) # Výpočet mediánu kapacit po 5 cyklech dle výrobců tapply(data5S$kap5, data5S$vyrobce, quantile, probs=0.5) # Obdobně lze určit další charakteristiky. # Pozor! Funkce pro výpočet šikmosti (skewness) a špičatosti (kurtosis) nejsou součástí základního R, najdete je v balíčku moments install.packages("moments") library(moments) # další charakteristiky -> var(), sd(), min(), max(), skewness(), kurtosis() ## Krabicový graf - angl. boxplot ##################################################### boxplot(b5) # Úprava grafu boxplot(b5, main="Kapacita po 5 cyklech (mAh)", xlab="Výrobce B", ylab="kapacita (mAh)") # Zabarvení boxplotu, zobrazení průměru boxplot(b5, main="Kapacita po 5 cyklech (mAh)", xlab="Výrobce B", ylab="kapacita (mAh)", col="grey") points(1, mean(b5,na.rm=TRUE), pch=3) # Horizontální boxplot boxplot(b5, main="Kapacita po 5 cyklech (mAh), výrobce B", horizontal=TRUE, xlab="kapacita (mAh)") # Další parametry boxplot(b5, main="Kapacita po 5 cyklech (mAh)", xlab="Výrobce B", ylab="kapacita (mAh)", col="grey", outline=FALSE, # nezobrazí odlehlá pozorování boxwex=0.5) # změní šířku krabice na 1/2 # Přidání přímky y=mean(b5) abline(h=mean(b5,na.rm=TRUE), col="red",lty=2) # Funkci boxplot lze využít nejen k vykreslení grafu boxplot(b5,plot=FALSE) ## Vícenásobný krabicový graf - angl. multiple boxplot ################################ boxplot(data5) boxplot(data5, boxwex=0.5, col="grey") boxplot(data5, boxwex=c(0.5,0.5,1.5,1.5), col=c("red","yellow","grey","blue")) boxplot(data5, boxwex=c(0.5,0.5,1.5,1.5), col=terrain.colors(4)) abline(h=1900,col="red",lty=2) # Mezera v boxplotu - myšlenka - třetí pozice bude vynechána boxplot(data5, at=c(1,2,4,5)) # Boxplot konstruován z dat ve standartním datovém formátu - s využitím funkce split pom=split(data5S$kap5, data5S$vyrobce) #do proměnné pom jsou uložena data ve formátu list, tj. seznam proměnných boxplot(pom) # Pořadí krabic dle přání s využitím funkce list pom=list(a5, b5, c5, d5) boxplot(pom) pom=list(d5, c5, b5, a5) boxplot(pom) ## Histogram ########################################################################## hist(a5) hist(a5,breaks=10) # Co dělají různé hodnoty parametru breaks s grafem? hist(a5, main="Histogram pro kapacitu akumulátorů po 5 cyklech, výrobce A", xlab="kapacita (mAh)", ylab="f(x)", col="blue", border="grey", labels=TRUE) # přidá absolutní četnosti daných kategorií ve formě popisků hist(a5, main="Histogram pro kapacitu akumulátorů po 5 cyklech, výrobce A", xlab="kapacita (mAh)", ylab="f(x)", col="blue", border="grey", labels=TRUE, # přidá absolutní četnosti daných kategorií ve formě popisků freq=FALSE) # změna měřítka na ose y --> f(x) lines(density(a5)) # připojí graf odhadu hustoty pravděpodobnosti # Generování hustoty normálního rozdělení xfit=seq(min(a5), max(a5), length=40) yfit=dnorm(xfit, mean=mean(a5), sd=sd(a5)) lines(xfit, yfit, col="black", lwd=2) ## QQ-graf - aneb grafický nástroj pro posouzení normality ############################ qqnorm(a5) qqline(a5) ####################################################################################### ## Více grafů do jednoho obrázku -> funkce layout nebo par ############################ ## Podívejte se na možnosti kombinování grafů - http://www.statmethods.net/advgraphs/layout.html ## Kombinace histogramu a boxplotu #################################################### pom=layout(mat = matrix(1:8,2,4, byrow=FALSE), height = c(2.5,1)) layout.show(pom) par(oma=c(2,2,3,2),mar=c(2,2,3,2)) hist(a5, main="Výrobce A", xlab="kapacita (mAh) po 5 cyklech", ylab="četnost", ylim=c(0,32), xlim=c(1730,2040)) boxplot(a5, horizontal=TRUE, ylim=c(1700,2040), boxwex=1.5) hist(b5, main="Výrobce B", xlab="kapacita (mAh) po 5 cyklech", ylab="četnost", ylim=c(0,32), xlim=c(1730,2040)) boxplot(b5, horizontal=TRUE, ylim=c(1700,2040), boxwex=1.5) hist(c5, main="Výrobce C", xlab="kapacita (mAh) po 5 cyklech", ylab="četnost", ylim=c(0,32), xlim=c(1730,2040)) boxplot(c5, horizontal=TRUE, ylim=c(1700,2040), boxwex=1.5) hist(d5, main="Výrobce D", xlab="kapacita (mAh) po 5 cyklech", ylab="četnost", ylim=c(0,32), xlim=c(1730,2040)) boxplot(d5, horizontal=TRUE, ylim=c(1700,2040), boxwex=1.5) mtext("Histogramy a boxploty pro jednotlivé výrobce po 5 cyklech", cex = 1.5, outer=TRUE, side=3) # Pro pokročilé - pomocí for-cyklu pom=layout(mat = matrix(1:8,2,4, byrow=FALSE), height = c(2.5,1)) layout.show(pom) par(oma=c(2,2,3,2), mar=c(2,2,3,2)) for (i in 1:4){ hist(data5[,i], main=paste("Výrobce",colnames(data5)[i]), xlab="kapacita (mAh) po 5 cyklech", ylab="četnost", xlim=c(min(data5,na.rm=TRUE), max(data5,na.rm=TRUE)), ylim=c(0,32)) boxplot(data5[,i], horizontal=TRUE, ylim=c(min(data5,na.rm=TRUE), max(data5,na.rm=TRUE)), boxwex=1.5) } mtext("Histogramy a boxploty pro jednotlivé výrobce po 5 cyklech", cex = 1.5, outer=TRUE, side=3) ## Kombinace histogramu a QQ-plotu #################################################### pom=layout(mat = matrix(1:8,2,4, byrow=FALSE), height = c(2,1.5)) layout.show(pom) par(oma=c(2,2,3,2), mar=c(2,2,3,2)) for (i in 1:4){ hist(data5[,i], main=paste("Výrobce",colnames(data5)[i]), xlab="kapacita (mAh) po 5 cyklech", ylab="četnost", xlim=c(min(data5,na.rm=TRUE), max(data5,na.rm=TRUE)), ylim=c(0,0.037), freq=FALSE) lines(density(data5[,i], na.rm=TRUE)) xfit=seq(min(data5[,i], na.rm=TRUE), max(data5[,i], na.rm=TRUE), length=40) yfit=dnorm(xfit, mean=mean(data5[,i], na.rm=TRUE), sd=sd(data5[,i], na.rm=TRUE)) lines(xfit, yfit, col="blue", lty=2) qqnorm(data5[,i]) qqline(data5[,i]) } mtext("Histogramy a QQ-ploty pro jednotlivé výrobce po 5 cyklech", cex = 1.5, outer=TRUE, side=3) ####################################################################################### ## 7. Identifikace odlehlých pozorování (a jejich odstranění z dat. rámce)############# ## S individuálním posouzením, jak naložit s odlehlými hodnotami ###################### # Seřázení datového rámce aku.st podle proměnných aku.st$vyrobce a data.st$kap5 data5S=data5S[with(data5S, order(data5S$vyrobce,data5S$kap5)),] # V RkWardu lze pro seřazení dat použít záložku Data / Sort data # Zápis proměnné do nového sloupce, v němž budeme odstraňovat vybraná odlehlá pozorování data5S$kap5.bez=data5S$kap5 # Následně lze pohodlně "ručně" odstranit (popř. opravit) vybraná odlehlá pozorování # Sofistikovanější metody -> např. vnitřní hradby, mnohorozměrná detekce odlehlých hodnot...