Data mining på svenska och maskininlärnings imputation
Jako Rostami
2020-05-20
Introduktion
Välkommen till denna genomgång av hur man minar, extrahearar, transformerar, visualiserar och imputerar data.
Datasetet kommer från https://openmv.net/info/travel-times och jag rekommenderar att du klickar där för att gå igenom beskrivningen av datasetet och variablerna innan du påbörjar detta. Det är en kort sammanfattning bara. Denna kod finns även på min github https://github.com/jakorostami
Du kommer att få lära dig nedan:
- Skriva funktioner med strängar som argument
- Hämta data och upptäcka mönster
- Enkel kod med dplyr och tidyverse
- Visualisering med ggplot
- Transformera och manipulera data
- Extrahera ny information från existerande
- Imputationsmetoder via regression och random forest såväl som logisk
- Skriva for loops
- Skriva nested ifelse funktioner
- Korrelationsanalys
- Hur man kodar =)
Då kör vi igång!
Datainhämtning
#Hämta hem alla paket
library(gridExtra)
library(corrplot)
library(scales)
library(RColorBrewer)
library(tidyverse)
library(ggplot2)
library(dplyr)
library(stats)
library(viridis)
library(hrbrthemes)
library(tidyr)
library(factoextra)
library(readxl)
library(knitr)
library(lubridate)
library(plyr)
#Vi börjar med att definiera en funktion för att söka efter NAs
na_list <- function(x){
if (sum(is.na(x)) > 0) {
sprintf("Du har %1.0f NAs. Skriv summary() för att se vilka kolumner de är i",
sum(is.na(x)))
} else {
print("Inga värden saknas =)")
}
}Vi hämtar hem data
##Hämta data
df <- read.csv2("travel-times.csv", sep=",", dec=".")
##Se i vilket format variablerna är i
str(df) #Använd glimpse() om du vill vara cool =)## 'data.frame': 205 obs. of 13 variables:
## $ Date : Factor w/ 111 levels "1/2/2012","1/3/2012",..: 4 4 3 3 2 2 1 1 53 52 ...
## $ StartTime : Factor w/ 123 levels "06:04","06:21",..: 75 48 70 29 122 31 106 16 33 96 ...
## $ DayOfWeek : Factor w/ 5 levels "Friday","Monday",..: 1 1 5 5 4 4 2 2 1 3 ...
## $ GoingTo : Factor w/ 2 levels "GSK","Home": 2 1 2 1 2 1 2 1 1 2 ...
## $ Distance : num 51.3 51.6 51.3 49.2 51.1 ...
## $ MaxSpeed : num 127 130 127 132 136 ...
## $ AvgSpeed : num 78.3 81.8 82 74.2 83.4 84.5 82.9 77.5 80.9 70.6 ...
## $ AvgMovingSpeed: num 84.8 88.9 85.8 82.9 88.1 88.8 87.3 85.9 88.3 78.1 ...
## $ FuelEconomy : Factor w/ 25 levels "","-","10.05",..: 1 1 1 1 1 1 2 2 17 17 ...
## $ TotalTime : num 39.3 37.9 37.5 39.8 36.8 36.8 37.2 37.9 39.3 43.5 ...
## $ MovingTime : num 36.3 34.9 35.9 35.6 34.8 35 35.3 34.3 36 39.3 ...
## $ Take407All : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ Comments : Factor w/ 23 levels "","Accident at 403/highway 6; detour along Dundas",..: 1 1 1 1 1 1 1 1 1 1 ...#Har vi NAs?
na_list(df) #Skall säga att inga värden saknas## [1] "Inga värden saknas =)"Kolla i datasetet så ser du att saknade värden inte representeras av NA
Vi hämtar hem datan på nytt genom att definiera NA för saknade värden
df <- read.csv2("travel-times.csv", sep=",", dec=".", na.strings = c("","-"))
#Har vi NAs?
na_list(df)## [1] "Du har 200 NAs. Skriv summary() för att se vilka kolumner de är i"Jämför med en påhittad matris utan missing data
x <- matrix(rnorm(16*21,0,1),16,21)
na_list(x)## [1] "Inga värden saknas =)"#Kolla avvikande outliers och se var NAs befinner sig
summary(df)## Date StartTime DayOfWeek GoingTo Distance
## 1/2/2012 : 2 08:11 : 5 Friday :27 GSK :105 Min. :48.32
## 1/3/2012 : 2 08:22 : 5 Monday :39 Home:100 1st Qu.:50.65
## 1/4/2012 : 2 16:15 : 5 Thursday :44 Median :51.14
## 1/6/2012 : 2 17:17 : 5 Tuesday :48 Mean :50.98
## 10/11/2011: 2 17:24 : 5 Wednesday:47 3rd Qu.:51.63
## 10/12/2011: 2 07:19 : 4 Max. :60.32
## (Other) :193 (Other):176
## MaxSpeed AvgSpeed AvgMovingSpeed FuelEconomy
## Min. :112.2 Min. : 38.10 Min. : 50.30 Min. : 7.810
## 1st Qu.:124.9 1st Qu.: 68.90 1st Qu.: 76.60 1st Qu.: 8.370
## Median :127.4 Median : 73.60 Median : 81.40 Median : 8.520
## Mean :127.6 Mean : 74.48 Mean : 81.98 Mean : 8.691
## 3rd Qu.:129.8 3rd Qu.: 79.90 3rd Qu.: 86.00 3rd Qu.: 8.970
## Max. :140.9 Max. :107.70 Max. :112.10 Max. :10.050
## NA's :19
## TotalTime MovingTime Take407All
## Min. :28.2 Min. :27.10 No :170
## 1st Qu.:38.4 1st Qu.:35.70 Yes: 35
## Median :41.3 Median :37.60
## Mean :41.9 Mean :37.87
## 3rd Qu.:44.4 3rd Qu.:39.90
## Max. :82.3 Max. :62.40
##
## Comments
## Backed up at Bronte : 2
## Rain, rain, rain : 2
## Accident at 403/highway 6; detour along Dundas: 1
## Accident blocked 407 exit : 1
## Accident: backup from Hamilton to 407 ramp : 1
## (Other) : 17
## NA's :181Vi ser att några kolumner har outliers. Det är de som befinner sig utanför 3e kvartilerna.
Här kommer vi inte att använda boxplots eftersom att en sannolikhetsfördelning då blir dold.
Histogram kommer att användas med tätheter för varje histogram.
Visuell analys
p <- df %>%
ggplot(aes(Distance, color=Take407All, fill=Take407All)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
facet_grid(Take407All ~ .)
#Här ska du fråga dig själv. Är det en outlier som skall bort eller har den ett praktiskt värde?
#Testa att skriva summary(df$Distance) och sd(df$Distance)
#Utgår man från ett konfidensintervall så är det extremt osannolikt att en sådan observation inträffar
#Det är en sannolikhet om ca 0.0000000000785% att observera ett avstånd om 60.32 km
#I sådana fall är det bra att kolla vilket datum det inträffade och kolla upp nyheter
#om trafiken eller förhållanden på vägen. Dock så är detta dataset för 1 person endast.
#Jag lämnar kvar den som ett eget beslut du får ta. Hint: Är den representabel för bilturerna?
p## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
p1 <- df %>%
ggplot(aes(MaxSpeed, color=Take407All, fill=Take407All)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
facet_grid(Take407All ~ .)
p1## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
p2 <- df %>%
ggplot(aes(AvgSpeed, color=Take407All, fill=Take407All)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
facet_grid(Take407All ~ .)
p2## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Vi kommer att jämföra denna senare
ja <- subset(df, df$Take407All == "Yes")
nej <- subset(df, df$Take407All == "No")
fuelmean <- c(mean(nej$FuelEconomy, na.rm=TRUE),
mean(ja$FuelEconomy, na.rm=TRUE))
fuelmean <- as.data.frame(fuelmean)
fuelmean$Take407All <- c("No", "Yes")
colnames(fuelmean) <- c("Mean", "Take407All")
fuelmean## Mean Take407All
## 1 8.764481 No
## 2 8.335000 Yesp3 <- df %>%
ggplot(aes(FuelEconomy, color=Take407All, fill=Take407All)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
geom_vline(data=fuelmean, aes(xintercept = Mean), linetype = "dashed") +
facet_grid(Take407All ~ .)
p3## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 19 rows containing non-finite values (stat_bin).## Warning: Removed 19 rows containing non-finite values (stat_density).
#Sätt ihop allt till ett fönster (du bör zooma/maximera fönstret)
grid.arrange(p,p1,p2,p3, ncol=4)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 19 rows containing non-finite values (stat_bin).
## Warning: Removed 19 rows containing non-finite values (stat_density).
#Låt oss se om vi kan studera variablerna TotalTime och MovingTime.
head(df)## Date StartTime DayOfWeek GoingTo Distance MaxSpeed AvgSpeed
## 1 1/6/2012 16:37 Friday Home 51.29 127.4 78.3
## 2 1/6/2012 08:20 Friday GSK 51.63 130.3 81.8
## 3 1/4/2012 16:17 Wednesday Home 51.27 127.4 82.0
## 4 1/4/2012 07:53 Wednesday GSK 49.17 132.3 74.2
## 5 1/3/2012 18:57 Tuesday Home 51.15 136.2 83.4
## 6 1/3/2012 07:57 Tuesday GSK 51.80 135.8 84.5
## AvgMovingSpeed FuelEconomy TotalTime MovingTime Take407All Comments
## 1 84.8 NA 39.3 36.3 No <NA>
## 2 88.9 NA 37.9 34.9 No <NA>
## 3 85.8 NA 37.5 35.9 No <NA>
## 4 82.9 NA 39.8 35.6 No <NA>
## 5 88.1 NA 36.8 34.8 No <NA>
## 6 88.8 NA 36.8 35.0 No <NA>#Vad kan de ha för syfte?
p4 <- df %>%
ggplot(aes(TotalTime, MovingTime)) +
geom_point(size=4) +
stat_smooth(aes(TotalTime, MovingTime), method="lm")
p4## `geom_smooth()` using formula 'y ~ x'
De verkar vara ko-linjära. Dvs att när den ena ökar så ökar den andra samt att båda beskriver tid men något dold verkar finnas mellan de.
TotalTime är total längd i minutes av bilresan och MovingTime är tiden bilen anses vara i rörelse.
Dvs tiden den kör exkl. trafik, olyckor eller när bilen står still.
Vi kan då få ut en ny variabel då bilen inte är i rörelse. Dvs nån slags väntetid.
Transformering
#Ta ut skillnaden mellan tidsvariablerna
res <- df$TotalTime - df$MovingTime
#Plotta skillnaderna med snitt och få fram deras egenskaper
plot(res, type="l", col="blue")
abline(h=mean(res), col="red")
res_eg <- data.frame(mean(res), sd(res), median(res))
colnames(res_eg) <- c("Medelvärde", "Standardavv.", "Median")
res_eg## Medelvärde Standardavv. Median
## 1 4.032683 3.048731 3.6Detta är en del av datawrangling som ett steg i datamining processen.
Vi lägger in den nya variabeln i vårt dataset.
#Skapa ny variabel utifrån två existerande
df <- df %>%
mutate(DelayTime = TotalTime - MovingTime)
#Vi kontrollerar fördelningen
p5 <- df %>%
ggplot(aes(DelayTime)) +
geom_histogram(aes(y=..density..),fill="navyblue") +
geom_density(alpha=.4, fill="black")
p5## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
##Här gäller samma som föregående gällande outliers. Låter den vara kvar som egen övning.
#Hint: Datatransformation - kvadratrot, normalisering, logaritmering etc.Låt oss kontrollera mot om man har tagit motorvägen eller inte.
p6 <- df %>%
ggplot(aes(DelayTime, color=Take407All, fill=Take407All)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
facet_grid(Take407All ~ .)
p6## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Vad är skillnaden mot om vi inte hade skapat en ny variabel?
j <- df %>%
ggplot(aes(TotalTime, color=Take407All, fill=Take407All)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
facet_grid(Take407All ~ .)
j1 <- df %>%
ggplot(aes(MovingTime, color=Take407All, fill=Take407All)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
facet_grid(Take407All ~ .)
grid.arrange(j,j1,p6, ncol=3)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Vi ser då tvetydiga grafer. Total tid minskar om man tar motorvägen men körtiden minskar samtidigt också. Ett enkelt svar borde vara att hastigheten ökar. Vi kan därför jämföra vår nya variabel mot AvgSpeed istället för att kolla på två ko-linjära variabler.
df %>%
ggplot(aes(AvgSpeed, DelayTime)) +
geom_point(size=4) +
stat_smooth(aes(AvgSpeed, DelayTime), method="loess")## `geom_smooth()` using formula 'y ~ x'
#Väntetiden sjunker ju högre snitthastighet man har! Visst är det självevident =)
#Men nu har du genom dataanalys framställt visualiseringar av dettaVi ska nu se vilka dagar som har högst median väntetid.
Ta median per veckodag och avrunda till 1 decimal.
medianer <- aggregate(DelayTime ~ DayOfWeek, df, mean)
medianer$DelayTime <- round(medianer$DelayTime, 1)
#Kör en boxplot med factor reverse
p7 <- df %>%
mutate(DayOfWeek = fct_reorder(DayOfWeek, DelayTime, .fun='median')) %>%
ggplot( aes(x=reorder(DayOfWeek, DelayTime), y=DelayTime, fill=DayOfWeek)) +
geom_boxplot() +
stat_summary(fun=median,
colour="yellow",
geom="point",
shape=18, size=3,
show.legend=FALSE) +
geom_text(data=medianer, aes(label=DelayTime, y= DelayTime - 0.6)) +
xlab("class") +
theme(legend.position="none") +
xlab("")
p7
Dock så är vi inte klara med vår missing data. Vår nyskapade variabel ska hjälpa oss med detta.
Först tar vi bort variabeln Comments eftersom att den inte bidrar med något av relevans för oss
df <- df[,-13]
#Kolla NAs - 19 NAs ska du ha
na_list(df)## [1] "Du har 19 NAs. Skriv summary() för att se vilka kolumner de är i"#Dessa skall finnas i FuelEconomy
summary(df)## Date StartTime DayOfWeek GoingTo Distance
## 1/2/2012 : 2 08:11 : 5 Friday :27 GSK :105 Min. :48.32
## 1/3/2012 : 2 08:22 : 5 Monday :39 Home:100 1st Qu.:50.65
## 1/4/2012 : 2 16:15 : 5 Thursday :44 Median :51.14
## 1/6/2012 : 2 17:17 : 5 Tuesday :48 Mean :50.98
## 10/11/2011: 2 17:24 : 5 Wednesday:47 3rd Qu.:51.63
## 10/12/2011: 2 07:19 : 4 Max. :60.32
## (Other) :193 (Other):176
## MaxSpeed AvgSpeed AvgMovingSpeed FuelEconomy
## Min. :112.2 Min. : 38.10 Min. : 50.30 Min. : 7.810
## 1st Qu.:124.9 1st Qu.: 68.90 1st Qu.: 76.60 1st Qu.: 8.370
## Median :127.4 Median : 73.60 Median : 81.40 Median : 8.520
## Mean :127.6 Mean : 74.48 Mean : 81.98 Mean : 8.691
## 3rd Qu.:129.8 3rd Qu.: 79.90 3rd Qu.: 86.00 3rd Qu.: 8.970
## Max. :140.9 Max. :107.70 Max. :112.10 Max. :10.050
## NA's :19
## TotalTime MovingTime Take407All DelayTime
## Min. :28.2 Min. :27.10 No :170 Min. : 0.000
## 1st Qu.:38.4 1st Qu.:35.70 Yes: 35 1st Qu.: 2.300
## Median :41.3 Median :37.60 Median : 3.600
## Mean :41.9 Mean :37.87 Mean : 4.033
## 3rd Qu.:44.4 3rd Qu.:39.90 3rd Qu.: 4.700
## Max. :82.3 Max. :62.40 Max. :25.800
## Nu till nästa steg. I datasetet har vi bara StartTime och inte EndTime. Vi har dock fyra viktiga variabler: StartTime, Distance, AvgSpeed och TotalTime.
#Kontrollera att TotalTime är korrekt genom nedan beräkning
rkoll <- rep(0, length(df$TotalTime))
koll <- rep(0, length(df$TotalTime))
for(i in 1:length(rkoll)){
#Detta ska resultera i samma värden som i kolumnen TotalTime
#Tid = Färdavstånd / Hastighet och sedan gånger 60 för att representera minuter
rkoll[i] <- (df[i,5] / df[i,7])*60
#Nedan skall resultera i 0 för varje rad då rkoll - TotalTime skall vara samma
#med hänsyn till avrundningsfel
koll[i] <- round((rkoll[i] - df[i,10]), 1)
}
koll #2 rader återger inte 0 vi kommer till dessa senare## [1] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.1 0.0 0.1 0.0 0.0 0.0 0.0
## [16] 0.0 0.0 0.0 0.0 -0.1 0.0 0.0 0.0 0.0 0.0 -0.1 0.0 0.0 0.0 0.0
## [31] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.1 0.0 0.0 -0.1 0.0 0.0
## [46] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 -0.1 0.0
## [61] 0.0 0.1 -0.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## [76] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 -0.1 0.0 -0.1 0.0 0.1 0.0
## [91] 0.0 0.0 0.0 -0.1 0.0 0.0 -0.1 0.0 48.1 0.0 0.0 0.0 0.0 0.0 0.0
## [106] 0.0 0.0 0.0 0.0 -0.1 0.0 0.0 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## [121] 0.0 0.1 0.0 0.0 0.0 0.0 0.0 0.5 0.1 0.1 0.0 0.0 0.0 0.0 0.0
## [136] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.1 0.0 -0.1 0.0 0.0 -0.1
## [151] 0.0 0.1 0.0 -0.2 0.0 -0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## [166] 0.0 -0.1 0.1 0.0 0.0 0.0 0.1 0.0 -0.1 0.0 -0.1 0.0 0.0 0.0 0.0
## [181] 0.0 0.0 -0.1 0.0 0.0 0.0 0.0 0.0 0.0 -0.1 0.0 -7.6 0.0 0.0 0.0
## [196] 0.0 -0.1 -0.1 0.0 0.0 0.1 0.0 0.0 0.0 -0.1max(koll); min(koll)## [1] 48.1## [1] -7.6Vi ska sätta ihop Date och StartTime till en variabel strax.
#Formatera datumet till avläsbart format
df$Date <- strptime(as.character(df$Date), "%m/%d/%Y") %>%
format("%Y-%m-%d")
#Sätt ihop datum och starttid
df$DateTime <- as.POSIXct(paste(df$Date, df$StartTime), format="%Y-%m-%d %H:%M")
#Skapa en funktion för skapa en ny variabel EndTime
endtime <- function(z){
t <- z * 3600
return(t)
}
#Transformera TotalTime till timmar
t2 <- rep(0, length(df$TotalTime))
for(i in 1:length(t2)){
t2[i] <- df[i,10] / 60
}
#Skapa nya variabeln EndDateTime och ta bort Date och StartTime - nya variablerna används senare
df$EndTime <- df[,14] + endtime(t2)
df <- df[,-c(1:2)]
#Hur stor del av datasetet är NA? Svar: 0.7%
dim(df)## [1] 205 13sprintf("%3.1f procent av datasetet är NA", (sum(is.na(df))/(205*13))*100)## [1] "0.7 procent av datasetet är NA"Imputation
För att använda imputationstekniker måste vi studera korrelation mellan variabler.
Att studera korrelationer mellan mellan variabler kan ge oss information om hur dessa står i relation till varandra.
Imputation är en teknik där man uppskattar det/de saknade värdena i ett dataset baserat på de faktiska värdena i andra variabler och/eller observationer i datasetet.
Du har säkert utfört en imputationsmetod själv - en typ som går under logisk imputation.
Säg att din hyresgäst har knappat in att han/hon har 33 barn i sin ansökan.
Med logisk imputation tar du bort en trea och låter det vara 3 barn.
#Skapa en korrelationsmatris
korr <- df[,c(3:9,11)] %>%
as.matrix() %>%
cor(use="pairwise.complete.obs")
#Plotta korrelationerna
#corrplot(korr, method="shade", order="hclust", col=brewer.pal(n=9, name="YlOrRd"))
corrplot(korr, method="number",
order="hclust",
col=brewer.pal(n=9, name="PuOr"),
type="upper",
number.cex=.6,
tl.col="black",
tl.cex=.8)
Vi verkar inte kunna använda oss av någon imputationsmetod sett till FuelEconomy.
Den har svaga korrelationer och FuelEconomy verkar påverkas av dolda faktorer.
Vi kan för genomgångens skull med ett diagnostest se om vi kan tillämpa en imputationsmetod.
1. Medelstort dataset (inom statistisk analys är ca 200 observationer medelstort)
2. Svaga korrelationer mellan variabler
3. ca 9% data saknas i variabeln FuelEconomy och 0.7% av hela datasetet
4. Våra saknade värden, genom att kolla på datasetet, verkar inte saknas efter någon slumpmässig följd. Data saknas i följd om man kollar på datumen.
Dvs det är olika dagar och olika datum men de sker i följd.
Punkt 1-4 talar för att vi kan ta bort dessa värden helt men vi kör en regressionsimputation.
avgspeedmean <- aggregate(AvgSpeed ~ DayOfWeek, df, mean)
movingtmean <- aggregate(MovingTime ~ DayOfWeek, df, mean)
modfuel <- lm(FuelEconomy ~ AvgSpeed + DayOfWeek + MovingTime - 1, df)
summary(modfuel)##
## Call:
## lm(formula = FuelEconomy ~ AvgSpeed + DayOfWeek + MovingTime -
## 1, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0589 -0.3052 -0.1151 0.3144 1.2067
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## AvgSpeed 0.01602 0.00574 2.791 0.00582 **
## DayOfWeekFriday 5.12079 0.91246 5.612 7.50e-08 ***
## DayOfWeekMonday 5.09848 0.90712 5.621 7.19e-08 ***
## DayOfWeekThursday 5.15739 0.90016 5.729 4.20e-08 ***
## DayOfWeekTuesday 4.98699 0.91065 5.476 1.45e-07 ***
## DayOfWeekWednesday 5.03186 0.91235 5.515 1.20e-07 ***
## MovingTime 0.06391 0.01359 4.703 5.10e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4771 on 179 degrees of freedom
## (19 observations deleted due to missingness)
## Multiple R-squared: 0.9971, Adjusted R-squared: 0.997
## F-statistic: 8820 on 7 and 179 DF, p-value: < 2.2e-16avgspeedmean## DayOfWeek AvgSpeed
## 1 Friday 81.65926
## 2 Monday 73.19744
## 3 Thursday 74.36591
## 4 Tuesday 73.78125
## 5 Wednesday 72.22979movingtmean## DayOfWeek MovingTime
## 1 Friday 35.11481
## 2 Monday 38.14615
## 3 Thursday 37.41818
## 4 Tuesday 38.42708
## 5 Wednesday 39.08511avgmovspeed <- data.frame(avgspeedmean, movingtmean)
avgmovspeed <- avgmovspeed[,-3]
avgmovspeed## DayOfWeek AvgSpeed MovingTime
## 1 Friday 81.65926 35.11481
## 2 Monday 73.19744 38.14615
## 3 Thursday 74.36591 37.41818
## 4 Tuesday 73.78125 38.42708
## 5 Wednesday 72.22979 39.08511idx <- which(is.na(df$FuelEconomy))
nalistan <- df[idx,c(1,5,9,7)]
#Fr 81.65926, Mon 73.19744, Thu 74.36591, Tue 73.78125, Wed 72.22979 avgspeed
#Fr 35.11481, Mon 38.14615, Thu 37.41818, Tue 38.42708, Wed 39.08511 movingtime
#Skapa en funktion med regressionekvationen
fuelreg <- function(x,y,z){
u = z ##Om du inte vill använda for-loop i nästa steg använd u <- deparse(substitute(z))
#Vi kör en nested ifelse
d = ifelse(u == "Friday", 5.12079,
ifelse(u == "Monday", 5.09848,
ifelse(u == "Thursday", 5.15739,
ifelse(u == "Tuesday", 4.98699,
ifelse(u == "Wednesday", 5.03186, 0)
)
)
)
)
#Och lägger in regressionsekvationen
t <- 0.01602*x + 0.06391*y + 1*d
return(t)
}
#Skapa en tom vektor
fu <- rep(0, 19)
#For loop för att köra regressionsimputation
for(i in 1:19){
fu[i] <- fuelreg(nalistan[i,2], nalistan[i,3], nalistan[i,1])
}
#Ersätt NAs i subsettet
nalistan$FuelEconomy <- round(fu, 2)
#Ersätt NAs i datasetet
df[idx,7] <- nalistan$FuelEconomy
##Skapa en egen dataframe för FuelEconomy som skall jämföras med random forest imputation senare
fu_eco <- df$FuelEconomyNu kan vi kolla skillnaden på före och efter imputationen.
impja <- subset(df, df$Take407All == "Yes")
impnej <- subset(df, df$Take407All == "No")
fuelmean2 <- c(mean(impnej$FuelEconomy), mean(impja$FuelEconomy))
fuelmean2 <- data.frame(fuelmean2)
fuelmean2$Take407All <- c("No", "Yes")
colnames(fuelmean2) <- c("Mean", "Take407All")
fuelmean2## Mean Take407All
## 1 8.758 No
## 2 8.338 Yesfuelmean2 <- rbind(fuelmean2,fuelmean)
tt <- c("Efter", "Efter", "Före","Före")
fuelmean2$Tid <- tt
colnames(fuelmean2) <- c("Mean", "Take407All", "Tid")
fuelmean2## Mean Take407All Tid
## 1 8.758000 No Efter
## 2 8.338000 Yes Efter
## 3 8.764481 No Före
## 4 8.335000 Yes Förep3## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 19 rows containing non-finite values (stat_bin).## Warning: Removed 19 rows containing non-finite values (stat_density).
px <- df %>%
ggplot(aes(FuelEconomy, color=Take407All, fill=Take407All)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
geom_vline(data=fuelmean2, aes(xintercept = Mean,
linetype=Tid)) +
facet_grid(Take407All ~ .)
px ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
grid.arrange(p3,px, ncol=2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 19 rows containing non-finite values (stat_bin).
## Warning: Removed 19 rows containing non-finite values (stat_density).## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Det gjorde knappt någon skillnad som vi också kan se från siffrorna
fuelmean2## Mean Take407All Tid
## 1 8.758000 No Efter
## 2 8.338000 Yes Efter
## 3 8.764481 No Före
## 4 8.335000 Yes FöreNu ska vi kontrollera tidskontrollen vi gjorde när vi kontrollerade TotalTime
Det är en logisk imputation via matematisk derivering
koll[99]## [1] 48.1df[99,]## DayOfWeek GoingTo Distance MaxSpeed AvgSpeed AvgMovingSpeed FuelEconomy
## 99 Thursday GSK 50.63 125.6 38.5 103.8 7.97
## TotalTime MovingTime Take407All DelayTime DateTime
## 99 30.8 29.3 Yes 1.5 2011-10-06 08:22:00
## EndTime
## 99 2011-10-06 08:52:48df[99,13] - df[99,12] #Time difference of 30.8 min## Time difference of 30.8 mins30.8/60 #tidkonvertering efter 60min = 1 tim## [1] 0.513333350.63/0.51333 #AvgSpeed 98.63 km/h ## [1] 98.6305150.63 / 98.63051 * 60 #30.8 avrundat - Korrekt!## [1] 30.7998#Detta innebär att vi har deriverat den verkliga snitthastigheten
#Och kan då imputera det felaktiga värdet utan någon imputationsmetod!
df[99,5] <- round(50.63 / (30.8/60), 1)
df[99,5]## [1] 98.6#Gör nu likadant med minimum värdet på -7.6 som ovan.Hints
Här visar jag små utdrag på vad man kan göra med outliers.
#I fallet med Distance
#Vi testar logaritmering
p <- df %>%
ggplot(aes(log(Distance))) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black")
p #Det hjälper inte## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Vi testar normalisering (standardisering av normalfördelning)
df$Distance <- scale(df$Distance)
p <- df %>%
ggplot(aes(Distance)) +
geom_histogram(aes(y=..density..), binwidth = 1) +
geom_density(alpha=.4, fill="black")
p #Det hjälpte inte heller
#Vi tar bort standardiseringen
df$Distance <- df$Distance * attr(df$Distance, 'scaled:scale') + attr(df$Distance, 'scaled:center')
which(df$Distance > 60) #Rad 45## [1] 45#Ta bort rad 45 och plotta igen
df <- df[-45,]
p <- df %>%
ggplot(aes(Distance)) +
geom_histogram(aes(y=..density..), binwidth = 1) +
geom_density(alpha=.4, fill="black")
p #Vi ser att vi har en bimodal fördelning efter flera transformationer
#Låt oss bekräfta detta genom att skala igen
df$Distance <- scale(df$Distance)
p <- df %>%
ggplot(aes(Distance)) +
geom_histogram(aes(y=..density..), binwidth = 1) +
geom_density(alpha=.4, fill="black")
p #Vi ser att vi har en bimodal fördelning efter flera transformationer
#Om vi plottar mot om man tar motorvägen eller inte för att upptäcka underliggande
p <- df %>%
ggplot(aes(Distance, color=Take407All, fill=Take407All)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
facet_grid(Take407All ~ .)
p #Vi ser att det är när man inte tar motorvägen som vi hittar två toppar## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
df$Distance <- df$Distance * attr(df$Distance, 'scaled:scale') + attr(df$Distance, 'scaled:center')
p <- df %>%
ggplot(aes(Distance, color=Take407All, fill=Take407All)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
facet_grid(Take407All ~ GoingTo)
p #Vi ser att efter splitten är det GoingTo GSK när man inte tar motorvägen som toppar bimodalt## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
bimodal <- df %>%
group_by(Take407All, GoingTo, Distance, DateTime, EndTime) %>%
filter(Take407All == "No", Distance < 50.5)
summary(bimodal)## DayOfWeek GoingTo Distance.V1 MaxSpeed AvgSpeed
## Friday : 6 GSK :38 Min. :48.32000 Min. :119.0 Min. :52.70
## Monday :10 Home: 1 1st Qu.:49.04000 1st Qu.:123.8 1st Qu.:65.65
## Thursday : 6 Median :49.11000 Median :124.8 Median :69.40
## Tuesday : 8 Mean :49.10692 Mean :125.8 Mean :68.44
## Wednesday: 9 3rd Qu.:49.16000 3rd Qu.:128.0 3rd Qu.:72.40
## Max. :50.18000 Max. :134.8 Max. :77.50
## AvgMovingSpeed FuelEconomy TotalTime MovingTime Take407All
## Min. :63.10 Min. :8.330 Min. :37.90 Min. :34.30 No :39
## 1st Qu.:73.95 1st Qu.:8.450 1st Qu.:40.70 1st Qu.:36.95 Yes: 0
## Median :77.70 Median :8.500 Median :42.40 Median :37.90
## Mean :76.43 Mean :8.579 Mean :43.43 Mean :38.76
## 3rd Qu.:79.70 3rd Qu.:8.540 3rd Qu.:44.90 3rd Qu.:39.90
## Max. :85.90 Max. :9.760 Max. :55.90 Max. :46.90
## DelayTime DateTime EndTime
## Min. :2.100 Min. :2011-07-25 08:06:00 Min. :2011-07-25 08:51:42
## 1st Qu.:3.100 1st Qu.:2011-08-09 20:14:00 1st Qu.:2011-08-09 20:58:33
## Median :4.500 Median :2011-08-24 07:59:00 Median :2011-08-24 08:49:18
## Mean :4.677 Mean :2011-09-13 13:13:53 Mean :2011-09-13 13:57:19
## 3rd Qu.:5.300 3rd Qu.:2011-09-16 07:12:30 3rd Qu.:2011-09-16 08:01:39
## Max. :9.300 Max. :2012-01-04 07:53:00 Max. :2012-01-04 08:32:48p <- bimodal %>%
ggplot(aes(Distance, color=GoingTo, fill=GoingTo)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
facet_grid(Take407All ~ DayOfWeek)
p #Här ser vi att vi kan utesluta GoingTo Home## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Groups with fewer than two data points have been dropped.
bimodal2 <- df %>%
group_by(Take407All, GoingTo, Distance, DateTime, EndTime) %>%
filter(Take407All == "No", Distance < 51, GoingTo == "GSK")
summary(bimodal2)## DayOfWeek GoingTo Distance.V1 MaxSpeed AvgSpeed
## Friday :6 GSK :38 Min. :48.32000 Min. :119.0 Min. :52.70
## Monday :9 Home: 0 1st Qu.:49.04000 1st Qu.:123.8 1st Qu.:65.62
## Thursday :6 Median :49.10500 Median :125.0 Median :69.25
## Tuesday :8 Mean :49.07868 Mean :125.9 Mean :68.38
## Wednesday:9 3rd Qu.:49.15750 3rd Qu.:128.2 3rd Qu.:72.60
## Max. :49.25000 Max. :134.8 Max. :77.50
## AvgMovingSpeed FuelEconomy TotalTime MovingTime Take407All
## Min. :63.10 Min. :8.330 Min. :37.90 Min. :34.30 No :38
## 1st Qu.:73.92 1st Qu.:8.450 1st Qu.:40.55 1st Qu.:36.92 Yes: 0
## Median :77.50 Median :8.500 Median :42.50 Median :38.00
## Mean :76.35 Mean :8.564 Mean :43.46 Mean :38.78
## 3rd Qu.:79.80 3rd Qu.:8.540 3rd Qu.:44.95 3rd Qu.:39.90
## Max. :85.90 Max. :9.760 Max. :55.90 Max. :46.90
## DelayTime DateTime EndTime
## Min. :2.100 Min. :2011-07-25 08:06:00 Min. :2011-07-25 08:51:42
## 1st Qu.:3.100 1st Qu.:2011-08-09 14:14:30 1st Qu.:2011-08-09 15:01:07
## Median :4.600 Median :2011-08-23 19:56:30 Median :2011-08-23 20:42:27
## Mean :4.682 Mean :2011-09-11 08:42:12 Mean :2011-09-11 09:25:40
## 3rd Qu.:5.300 3rd Qu.:2011-09-13 01:54:30 3rd Qu.:2011-09-13 02:47:37
## Max. :9.300 Max. :2012-01-04 07:53:00 Max. :2012-01-04 08:32:48#Vilka datum har bilturerna skett?
bimodal2 %>%
ggplot(aes(as.Date(DateTime))) +
geom_histogram() +
scale_x_date(labels = date_format("%Y-%m-%d"))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#I början på januari försvinner detta fenomen. Nu känner vi inte till personens
#livsföring men man kan anta att det har varit en väg som har varit förekommande
#under en tillfällig men återkommande period. Sedan har personen anpassat sina
#bilturer till GSK efter januari 2012 som motsvarar den nuvarande körningen.
#Vi tar bort de gamla turerna också
idx2 <- which(df$Distance < 50.5)
df2 <- df[-idx2,]
pdf2 <- df2 %>%
ggplot(aes(Distance, color=Take407All, fill=Take407All)) +
geom_histogram(aes(y=..density..)) +
geom_density(alpha=.4, fill="black") +
facet_grid(Take407All ~ .)
pdf2 ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#De bimodala topparna är nu borta med mindre datapunkter dock.Imputation med slumpmässiga beslutsträd kallat Random Forest
Det sista jag vill visa är hur man kan använda statistiska algoritmer (machine learning) för att imputera saknade värden. Ett av de är random forest - slumpmässiga beslutsträd.
Observera att det inte finns någon definitiv metod för imputation och att det är något som skall gås igenom grundligt!
– Vad är random forest och hur funkar det? –
Bootstrapping
Bootstrapping är att ta slumpmässiga stickprov ur vårt ordinarie dataset för att skapa ett dataset lika stort som vårt ordinarie dataset. Vi kan ta samma stickprov flera gånger.Beslutsträd
Skapa ett beslutsträd från ditt bootstrappade dataset och använd endast en del av av dina variabler vid varje steg (gren). Om du har 4 variabler, använd 2 (slumpmässigt valda) vid första grenen och 2 vid sista grenen.Skapa en slumpmässig skog (dålig översättning av random forest) genom att repetera stegen ovan. Om du repeterar stegen ovan flera hundra gånger blir det en fin liten skog =)
#library(randomForest) Ta bort # för att ladda paketen
#library(missForest)
#Vi börjar att hämta om datan så att vi har NAs
df <- read.csv2("travel-times.csv", sep=",", dec=".", na.strings = c("","-"))
#Ta bort Comments
df <- df[,-13]
na_list(df)## [1] "Du har 19 NAs. Skriv summary() för att se vilka kolumner de är i"#För att reproducera kommer vi att sätta ett seed
set.seed(1234)
#Vi skapar en ny dataframe med missForest paketet imputation
df_imp <- missForest(df[,3:12], #Vi skippar datum och tid variablerna
ntree = 100, #Vi odlar 100 träd i varje skog
maxiter = 6, #Vi sätter en gräns till max 6 skogar för att imputera
verbose = TRUE #Följ iterationerna och se estimerade fel
)## missForest iteration 1 in progress...done!
## estimated error(s): 0.005090842 0
## difference(s): 1.778259e-07 0
## time: 0.05 seconds
##
## missForest iteration 2 in progress...done!
## estimated error(s): 0.005056097 0
## difference(s): 1.847217e-08 0
## time: 0.03 seconds
##
## missForest iteration 3 in progress...done!
## estimated error(s): 0.005099556 0
## difference(s): 1.872613e-08 0
## time: 0.03 seconds#Nu har vi fått imputerade värden där det saknades i FuelEconomy
#Vi frågar om det finns några NAs
anyNA(df_imp$ximp$FuelEconomy) #False## [1] FALSE#Jämför med dina imputerade värden från regressionsekvationen
imp_diff <- data.frame(fu_eco, df_imp$ximp$FuelEconomy)
colnames(imp_diff) <- c("FuelEco RegImp", "FuelEco rfImp")
#Granska de NAs rader som imputerats
View(imp_diff[idx,])
head(imp_diff[idx,])## FuelEco RegImp FuelEco rfImp
## 1 8.70 8.6024
## 2 8.66 8.8866
## 3 8.64 8.8032
## 4 8.50 8.6871
## 5 8.55 8.4591
## 6 8.58 8.5845#Se differenserna
imp_diff[idx,1] - imp_diff[idx,2]## [1] 0.0976 -0.2266 -0.1632 -0.1871 0.0909 -0.0045 0.1856 0.0686 -0.3778
## [10] -0.1507 0.1847 -0.4165 0.0708 0.1905 -0.0094 -0.3564 0.0968 -0.0816
## [19] -0.3529Det var allt för denna gång. Välkommen tillbaka!! =)