Kaggle assignment
Giulia Gualdi
12 maggio 2016
library(ggplot2) #visualization
library(ggthemes)#visualization
library(rpart) # ml algorithm:trees
library(randomForest) # ml algorithm: random forests
#library(caret) ml algorithms
library(dplyr) #data manipulation
library(lubridate) #date manipulation
library(caret)#machinelearning
library(rattle) #tree visualizationRead in test & training data
#setwd('/Users/giulia/datasci_course_materials/Assignments/assignment6')
data = read.csv('train.csv',na.strings=c(""))
test = read.csv('test.csv',na.strings=c(""))and check consistency between columns of data and test datasets
names(data)## [1] "AnimalID" "Name" "DateTime" "OutcomeType"
## [5] "OutcomeSubtype" "AnimalType" "SexuponOutcome" "AgeuponOutcome"
## [9] "Breed" "Color"names(test)## [1] "ID" "Name" "DateTime" "AnimalType"
## [5] "SexuponOutcome" "AgeuponOutcome" "Breed" "Color"length(names(data)) == length(names(test))+1## [1] FALSEthe two datasets contain a different number of columns and some the ‘ID’ columns are named differently. First thing, we fix this inconsistency
names(data)[1]='ID'In addition, from visual inspection, we can also see that the test dataset does not contain the ‘OutcomeSubtype’ column which means that we won’t use it as a feature.
data$OutcomeSubtype<- NULLNow we can check for consistency in terms of variable types usin the str() function and we can fix the inconsistency in the ID columns
data$ID = as.integer(data$ID)Exploratory data analysis
Let’s fist start with an overview of the data
summary(data)## ID Name DateTime
## Min. : 1 Max : 136 2015-08-11 00:00:00: 19
## 1st Qu.: 6683 Bella : 135 2015-11-17 00:00:00: 17
## Median :13365 Charlie: 107 2015-07-02 00:00:00: 13
## Mean :13365 Daisy : 106 2015-04-02 00:00:00: 11
## 3rd Qu.:20047 Lucy : 94 2014-08-31 09:00:00: 10
## Max. :26729 (Other):18460 2014-08-26 09:00:00: 9
## NA's : 7691 (Other) :26650
## OutcomeType AnimalType SexuponOutcome AgeuponOutcome
## Adoption :10769 Cat:11134 Intact Female:3511 1 year : 3969
## Died : 197 Dog:15595 Intact Male :3525 2 years : 3742
## Euthanasia : 1555 Neutered Male:9779 2 months: 3397
## Return_to_owner: 4786 Spayed Female:8820 3 years : 1823
## Transfer : 9422 Unknown :1093 1 month : 1281
## NA's : 1 (Other) :12499
## NA's : 18
## Breed Color
## Domestic Shorthair Mix : 8810 Black/White : 2824
## Pit Bull Mix : 1906 Black : 2292
## Chihuahua Shorthair Mix : 1766 Brown Tabby : 1635
## Labrador Retriever Mix : 1363 Brown Tabby/White: 940
## Domestic Medium Hair Mix: 839 White : 931
## German Shepherd Mix : 575 Brown/White : 884
## (Other) :11470 (Other) :17223The data summary shows that approx 60% of data contained in the training set is about dogs, and, as a general remark, we can guess that on potentially relevant features (like dog breed) there will be some cleaning to do.
In order to go a little deeper, let’s see how the outcomes are distributed for cats and dogs
ggplot(data, aes(OutcomeType,col=AnimalType,fill=AnimalType)) + geom_bar(aes(y = (..count..)/ sapply(PANEL, FUN=function(x) sum(count[PANEL == x])))) +facet_wrap(~ AnimalType)+ labs(y = 'Frequency',
x = 'Outcome',
title = 'OutcomeType: Cats vs Dogs') +
theme(axis.text.x = element_text(face="bold", angle=45))
From the bar plot it seems that cats are more likely to be transferred than dogs, who are instead more likely to be returned to their owners. Luckily, in both cases, the dreadful options (death & euthanasia) are not that frequent (less than 10% of occurrences in both categories)
In the following I will first examine animal-related features (e.g. Age, breed, etc) and then the ‘environmental’ variable, i.e. date-time
Animal-related features
Age
Let’s begin by examining ‘AgeuponOutcome’. My blind guess is that the age of the animal might have a huge influence on the outcome as people are in general more prone to adopt kittens/puppies, that very old and very young animals, being on the weaker side, might be more exposed to death & euthanasia outcomes, whereas, perhaps, adult animals might be more likely to return to their owners. Let’s start to see how many levels this variable has
levels(data$AgeuponOutcome)## [1] "0 years" "1 day" "1 month" "1 week" "1 weeks"
## [6] "1 year" "10 months" "10 years" "11 months" "11 years"
## [11] "12 years" "13 years" "14 years" "15 years" "16 years"
## [16] "17 years" "18 years" "19 years" "2 days" "2 months"
## [21] "2 weeks" "2 years" "20 years" "3 days" "3 months"
## [26] "3 weeks" "3 years" "4 days" "4 months" "4 weeks"
## [31] "4 years" "5 days" "5 months" "5 weeks" "5 years"
## [36] "6 days" "6 months" "6 years" "7 months" "7 years"
## [41] "8 months" "8 years" "9 months" "9 years"and if there is any missing value both in the training set
length(which(is.na(data$AgeuponOutcome)))## [1] 18and in the test set
length(which(is.na(test$AgeuponOutcome)))## [1] 6Since we have missing data in the test set as well we need to fill in the missing values. Since I want to alter the data the least possible, I remove the 18 rows corresponding to the missing values in the train dataset
data=data[-which(is.na(data$AgeuponOutcome)),]while I will fill in the missing values in the dataset only once I will have figured out which is the most meaningful labelling for the age feature.
I now proceed to see how the data is distributed in terms of years by grouping all the dates in days/weeks/months under the ‘O years’ label and creating a new numeric column ‘AgeinYears’ based on which we will sort the dataset (good for plotting purposes..)
data$AgeinYears = data$AgeuponOutcome
b = data$AgeinYears
ind_train = union(union(grep('month',b),grep('day',b)),grep('week',b))
data$AgeinYears[ind_train] = '0 years'
data$AgeinYears= as.numeric(gsub('year|years',"",data$AgeinYears))
data= data[order(data$AgeinYears),]..finally we can see how the ages are distributed among cats & dogs along with the outcomes
ggplot(data, aes(AgeinYears
,col=OutcomeType,fill=OutcomeType)) + geom_bar(aes(y = (..count..)/ sapply(PANEL, FUN=function(x) sum(count[PANEL == x])))) +facet_wrap(~ AnimalType)+ labs(y = 'Frequency',
x = 'Age of Outcome (years)',title = 'Ageof Outcome and OutcomeType: Cats vs Dogs') +
theme(axis.text.x = element_text(face="bold", angle=45))
Interestingly, while more than 65% of cats are less than one year upon outcome, only almost 30% of dogs are, thus indicating a somewhat less skewed dog age distribution in our sample. Also, it seems that the decrease in adoption likelihood for cats older than 1 years is more pronounced. Death seems almost exclusively an issue for kittens, while euthanasia is more evenly spread along all ages and animal categories.
We could aggregate ages further in order to reduce complexity and increase later interpretability. Based on my personal experience I break ages down in four categories: junior (<1 year), young adults (1-5 yrs), adults (5-10), mature (>10)
data$AgeCat = ifelse(data$AgeinYears == 0, 'junior',
ifelse(data$AgeinYears >=1 & data$AgeinYears <= 5, 'young adult',
ifelse(data$AgeinYears > 5 & data$AgeinYears <= 10,'adult','mature')))
data$AgeCat = factor(data$AgeCat,levels = c('junior','young adult','adult','mature'))and re-check the outcome distribution normalized over those categories
bycat = group_by(data,AnimalType,AgeCat,OutcomeType)
sumc = summarise(bycat,count = n())
ggplot(sumc, aes(x = AgeCat, y = count, fill = OutcomeType)) +
geom_bar(stat = 'identity', position = 'fill', colour = 'black') +
facet_wrap(~AnimalType) +
coord_flip() +
labs(y = 'percentage of animals',
x = 'Age',
title = 'Outcomes by Age Category: Cats & Dogs') 
Since the aggregated plot is still very informative, we proceed to add the new columns to the test dataset as well
test$AgeinYears = test$AgeuponOutcome
c = test$AgeinYears
ind_test = union(union(grep('month',c),grep('day',c)),grep('week',c))
test$AgeinYears[ind_test] = '0 years'
test$AgeinYears= as.numeric(gsub('year|years',"",test$AgeinYears))
test$AgeCat = ifelse(test$AgeinYears == 0, 'junior',
ifelse(test$AgeinYears >=1 & test$AgeinYears <= 5, 'young adult',
ifelse(test$AgeinYears > 5 & test$AgeinYears <= 10,'adult','mature')))
test$AgeCat = factor(test$AgeCat,levels = c('junior','young adult','adult','mature'))
test$AgeinYears=NULL
data$AgeinYears=NULLSince for classification I will use the AgeCat feature, I will now train a decision tree in order to fill in the missing values in the test dataset. In principle I could use the full dataset (6 missing values excluded) to train the tree, however I prefer to be on the safe side and use only the training dataset
ageform = AgeCat ~ AnimalType
agemodel = rpart(ageform,method="class", data=data)
print(agemodel)## n= 26711
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 26711 14821 junior (0.44513496 0.43416570 0.09250870 0.02819063)
## 2) AnimalType=Cat 11117 3758 junior (0.66195916 0.26580912 0.05253216 0.01969956) *
## 3) AnimalType=Dog 15594 6952 young adult (0.29056047 0.55418751 0.12100808 0.03424394) *agepred = predict(agemodel, newdata = test[which(is.na(test$AgeCat)),])
ageprediction = (colnames(agepred)[max.col(agepred,ties.method="first")])
test$AgeCat[which(is.na(test$AgeCat))]=agepredictionBreed
Next, we proceed to examine the breed feature. As a preliminary step I check that there are no missing values
which(is.na(data$Breed))## integer(0)which(is.na(test$Breed))## integer(0)In my opinion, the breed should be especially relevant to dogs, mainly with respect to adoption and return to owner. In fact, many people who opt for dog adoption might be somewhat sensitive to some specific features of the breed (e.g. size, energy levels, sociability) to the extent that those can be informative of how well suited the pet could be to their lifestyle. This is less of a problem for cats as their behaviour is, except for few breeds, less-breed related and more dependent on the indidual animal and their size-variations are not as significant as in dogs.
In few words, when it comes to adoption, my feeling is that most of the times it is not so much the specific breed that counts but rather some related features, e.g. size, energy and sociability. For this reason I am confident that the 1380 different breeds listed in the dataset can be drastically reduced.
length(levels(data$Breed))## [1] 1380From visual inspection of the entries of the Breed column my understanding is that we have three type of entries: 1) Breed 2) Breed mix; 3) Breed A/Breed B;
the difference between 2 and 3 being that in the case of 2 it is not possible to identify which cross occurred. If my guess is correct then, at first approximation, entries of kind 2 and 3 are essentially the same thing and can be collapsed into 2, and we can add a column to keep track of crosses vs pure breeds
data$AggBreed=sapply(as.character(data$Breed),function(x) ifelse(grepl('/',x),paste(strsplit(x, split = '/')[[1]][1],'Mix',sep=' '),x))
data$IsPureBreed = ifelse(grepl('Mix',data$AggBreed),0,1)Now, by looking at the unique values of AggBreed, we see that the Domestic cat is considered as a pure breed, but my guess is that those who want pure breed cats definitely do not think about the domestic cat! In addition, I believe that the indication about the length of the hair of the domestic cat is a bit too much of a detail for our analysis. For this reason I collapse all the variations of Domestic cat into a single one and classify it as mixed as, according to Wikipedia, the ‘Domestic cat’ is a term used to indicate cats without a specific breed
data$IsPureBreed[grep('Domestic',data$AggBreed)]=0
data$AggBreed[grep('Domestic',data$AggBreed)] = 'Domestic'
length(data$AggBreed)## [1] 26711Now we are left with 350 different breeds, which is less than one third of the initial amount, but still way too many. For later choice of the optimal optimization set, I want to see how many dog vs cat breeds we have, and gauge to some extent the relevance of the breed in the two animal categories.
aggregate(data=data, AggBreed ~ AnimalType, function(x) length(unique(x)))## AnimalType AggBreed
## 1 Cat 42
## 2 Dog 308so we see that we have 42 breed labels for cats vs 308 for dogs. In order to reduce further, the best would be to aggregate dog breeds by their ‘official’ (i.e. assigned by a relevant Kennel club) group label, as it is in general an indicator of the relevant macro-features of the breeds, but for the moment I will stick with what I have and perhaps try to include this information in later versions of this script. In order to proceed and see if any further reduction is possible, I want to see if an appreciable different pattern emerges for pure-breed vs mix
bycat = group_by(data,AnimalType,IsPureBreed,OutcomeType)
sumc = summarise(bycat,count = n())
ggplot(sumc, aes(x = IsPureBreed, y = count, fill = OutcomeType)) +
geom_bar(stat = 'identity', position = 'fill', colour = 'black') +
facet_wrap(~AnimalType) +
coord_flip() +
labs(y = 'percentage of animals',
x = 'Breed purity',
title = 'Outcomes by Age Category: Cats & Dogs') 
From the plot it seems, quite interestingly!, that whereas breed purity is positively correlated with adoption for cats, the inverse is true for dogs. However, these might well be simple sample fluctuations and it might actually well be that breed purity has no significant impact. Indeed, my gut feeling is that breed purity is not really relevant when it comes to shelter animals and that, as I already mentioned, other breed-related features are the crucial ones.
However, since there are too many breeds and no meaningful way to boil this feature down to few relevant classes without the help of an external database, for the moment the only breed info I will use in the classification algorithm is the one about purity.
Finally, let’s add the new columns to the test set as well and remove AggBreed from both datasets once it has been used
test$AggBreed=sapply(as.character(test$Breed),function(x) ifelse(grepl('/',x),paste(strsplit(x, split = '/')[[1]][1],'Mix',sep=' '),x))
test$IsPureBreed = ifelse(grepl('Mix',test$AggBreed),0,1)
test$IsPureBreed[grep('Domestic',test$AggBreed)]=0
#test$AggBreed[grep('Domestic',test$AggBreed)] = 'Domestic'
test$AggBreed=NULL
data$AggBreed=NULLSexUponOutcome
Let’s now examine the SexuponOutcome feature. Let’s first see if there is any missing entry
which(is.na(data$SexuponOutcome))## [1] 24232which(is.na(test$SexuponOutcome))## integer(0)Since there ise one missing entry in the data set, it seems to me reasonble to classify it as unknown
data$SexuponOutcome[which(is.na(data$SexuponOutcome))]='Unknown'At this point we have five labels:
levels(data$SexuponOutcome)## [1] "Intact Female" "Intact Male" "Neutered Male" "Spayed Female"
## [5] "Unknown"#data$SexuponOutcome[which(data$SexuponOutcome == '')] = 'Unknown'I now plot the outcome by Sex to see if any pattern emerges and if the variable needs to be broken down further
bycat = group_by(data,AnimalType,SexuponOutcome,OutcomeType)
sumc = summarise(bycat,count = n())
ggplot(sumc, aes(x = SexuponOutcome, y = count, fill = OutcomeType)) +
geom_bar(stat = 'identity', position = 'fill', colour = 'black') +
facet_wrap(~AnimalType) +
coord_flip() +
labs(y = 'percentage of animals',
x = 'SexUponOutcome',
title = 'Outcomes by Sex: Cats & Dogs') 
From the chart we see that for both animals the adoptions are concentrated in the neutered/spayed areas, and that somehow the sex of the animal does not seem too relevant. In order to check this, I will create a variable, ‘operated’ which groups intact and unknowns together vs spayed/neutered
data$operated =
ifelse(data$SexuponOutcome == 'Neutered Male' |
data$SexuponOutcome == 'Spayed Female','yes','no/maybe')Next, I want to check how the ‘operated’ feature and the age feature are correlated
bycat = group_by(data,AnimalType,AgeCat,operated)
sumc = summarise(bycat,count = n())
ggplot(sumc, aes(x = AgeCat, y = count, fill = operated)) +
geom_bar(stat = 'identity', position = 'fill', colour = 'black') +
facet_wrap(~AnimalType) +
coord_flip() +
labs(y = 'percentage of animals',
x = 'Age Category',
title = 'percentage of neutered/spayed animals by age category:
Cats & Dogs') 
From the chart we can see that even if some pattern might emerge, this is definitely not enough to consider the two features almost as ‘synonims’ so we will keep them both.
We now proceed to add the ‘operated’ column to the test set as well
test$SexuponOutcome[which(test$SexuponOutcome == '')] = 'Unknown'
test$operated = ifelse(test$SexuponOutcome == 'Neutered Male' |
test$SexuponOutcome == 'Spayed Female','yes','no/maybe')
data$operated = as.factor(data$operated)
test$operated = as.factor(test$operated)Color
The color of an animal might contribute its visual appeal and hence increase its chances to be adopted even though I do not expect this variable to have any relevance with respect to any outcome other than this and return to owner. Even in the latter cases I am prone to believe that the color will not be that relevant.
Let’s check first how many colors we have animal category
length(levels(data$Color))## [1] 366367, way too many! and how many per animal type
aggregate(data=data, Color ~ AnimalType, function(x) length(unique(x)))## AnimalType Color
## 1 Cat 146
## 2 Dog 262Conclusion:, I will not include the color feature in my classification algorithm because I do not believe that this feature is extremely relevant nor informative so It is not worth the effort of trying to reduce the hundreds of labels to less than then (my personal threshold for useful labels)
Name
Finally, the last animal-related feature is ‘Name’. I personally believe that the fact that an animal has a name or not has no whatsoever relevance to the outcome, as many adopters would nevertheless give their pets another name and for the other categories I do not see how it could be relevant. I will therefore neither manipulate this variable nor include it in the classification algorithm.
Environmental features
DateTime
The date/time feature might contain relevant information to guess the outcome of an animal. I would expect some seasonality dependence not only because humans might be more prone to get a pet during periods where they do not plan to go on holiday, but also because some animal deseases might depend on the climate as perhaps those insect-related. I might also guess that the chances of return to owner and adoption peak in the week ends and those of transfer or euthanasia during week days, with death, obviously, being evenly spread. Same goes for the hour of the day, as death can occurr at any time of the day whereas outcomes which involve human intervention will tend to be more concentrated during daily hours.
Ok let’s get to work!
which(is.na(data$DateTime))## integer(0)which(is.na(test$DateTime))## integer(0)data$month=month(data$DateTime)
data$day=wday(data$DateTime)
data$year = year(data$DateTime)
data$hour = hour(data$DateTime)
winter = c(12,1,2)
spring = c(3,4,5)
summer = c(6,7,8)
autumn = c(9,10,11)
weekend = c(6,7)
working = c(1,4)
morning=c(8,9,10,11)
earlymor = c(5,6,7)
midday = c(12,13,14,15)
afternoon =c(16,17,18)
evening = c(19,20,21)
night = c(22,23,0,1,2,3,4)
data$season = ifelse(is.element(data$month,winter),'winter',
ifelse(is.element(data$month,spring),'spring',
ifelse(is.element(data$month,autumn),'autumn','summer')))
data$daytype = ifelse(is.element(data$day,weekend),'weekend',
ifelse(is.element(data$day,working),'working','Friday'))
data$time = ifelse(is.element(data$hour,earlymor),'early morning',
ifelse(is.element(data$hour,morning),'morning',
ifelse(is.element(data$hour,midday),'midday',
ifelse(is.element(data$hour,evening),'evening','night'))))
data$season = as.factor(data$season)
data$daytype = as.factor(data$daytype)
data$time = as.factor(data$time)where I kept Friday as a special daytype as it is often a shorter working day and usually the mood is already week-end oriented so I felt it should be treated differently with respect to, say, a Tuesday
bycat = group_by(data,AnimalType,season,OutcomeType)
sumc = summarise(bycat,count = n())
ggplot(sumc, aes(x = season, y = count, fill = OutcomeType)) +
geom_bar(stat = 'identity', position = 'fill', colour = 'black') +
facet_wrap(~AnimalType) +
coord_flip() +
labs(y = 'percentage of animals',
x = 'Season',
title = 'Outcomes by season') 
Interestingly, the data show no sensational seasonal pattern. Nevertheless I will still keep this feature as I feel overlooking it might be a mistake. Let’s check if we have better luck with weekly patterns
bycat = group_by(data,AnimalType,daytype,OutcomeType)
sumc = summarise(bycat,count = n())
ggplot(sumc, aes(x = daytype, y = count, fill = OutcomeType)) +
geom_bar(stat = 'identity', position = 'fill', colour = 'black') +
facet_wrap(~AnimalType) +
coord_flip() +
labs(y = 'percentage of animals',
x = 'Type of day',
title = 'Outcomes by weekdays')
Nope, from our data we do not see striking weekly patters. I will nevertheless keep this variable and include it in the test set because it might be helpful at later stages.
test$month=month(test$DateTime)
test$day=wday(test$DateTime)
test$year = year(test$DateTime)
test$hour = hour(test$DateTime)
test$season = ifelse(is.element(test$month,winter),'winter',
ifelse(is.element(test$month,spring),'spring',
ifelse(is.element(test$month,autumn),'autumn','summer')))
test$daytype = ifelse(is.element(test$day,weekend),'weekend',
ifelse(is.element(test$day,working),'working','Friday'))
test$time = ifelse(is.element(test$hour,earlymor),'early morning',
ifelse(is.element(test$hour,morning),'morning',
ifelse(is.element(test$hour,midday),'midday',
ifelse(is.element(test$hour,evening),'evening','night'))))
test$season = as.factor(test$season)
test$daytype = as.factor(test$daytype)
test$time = as.factor(test$time)Summary of the preliminary analysis
After our various manipulations we ended up summarizing information into the new columns AgeCat,IsPureBreed,operated,season,daytipe. Breed,Color,DateTime,OutcomeSubtype will surely not be included in the classification algorithm.
Classification
Now we are ready to proceed to the classification. First thing I want to do is to split the data into training and validation, in order to be able to compare various models and algorithms.
set.seed(111)
ind = createDataPartition(data$OutcomeType,p=0.75)[[1]]#must access the first (and only) element
#of the list to be used as an index array
train = data[ind,]
valid = data[-ind,]First I want to create a benchmark model where the probability of each outcome is assigned based onits frequency in the sample. I do not use the logloss function to score the results as in the valid set I do not have probabilities of outcomes but the actual outcomes so it seems to me more meaningful to simply count the number of entries which the model got right and divide by the total
freq = table(data$OutcomeType)[1:5]/nrow(data)
bench = as.data.frame(matrix(ncol=2,nrow=nrow(valid)))
names(bench)=c('ID','OutcomeType')
bench$ID = valid$ID
outcomes = levels(data$OutcomeType)
p= c(freq[[1]],freq[[2]],freq[[3]],freq[[4]],freq[[5]])
sampleDist = function(n) { sample(x = outcomes, n, replace = T, prob = p)}
bench$OutcomeType=sampleDist(nrow(bench))
score = sum(bench$OutcomeType==as.character(valid$OutcomeType))/nrow(valid)
score## [1] 0.3132115Ok, so we see that our benchmark model scored pretty badly with only 30% of correct predictions, so hopefully it should not be that hard to do better!
As a second step I want to classify the outcomes using a decision tree. Even though in the end I will also use a random forest, I want to check how the two compare.
formula = OutcomeType~ AnimalType + AgeCat + IsPureBreed + operated + season + daytype + time
modeldt <-rpart(formula, method="class", data=train)
#print(modeldt)
fancyRpartPlot(modeldt)
according to the decision tree the most important variables in determining the outcome type are whether the animal has been neutered/spayed, its age and then the animal type. Interestingly we see that, according to the decision tree, if an animal has not been operated then transfer is the only option.
Despite using all our variables we see that the decision tree manages to distinguish only 3 out of 5 outcomes, not detecting Death and Euthanasia.
Let’s see how this is translated in terms of accuracy
pred = predict(modeldt, newdata = valid)
prediction = (colnames(pred)[max.col(pred,ties.method="first")])
acc = sum(prediction== as.character(valid$OutcomeType))/length(prediction)
acc## [1] 0.6417016With the decision tree we only got 64% of instances right. Let’s see if using a random forest we can do better
model1rf <- randomForest(formula, data=train)
print(model1rf)##
## Call:
## randomForest(formula = formula, data = train)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 36.38%
## Confusion matrix:
## Adoption Died Euthanasia Return_to_owner Transfer
## Adoption 7146 0 0 448 483
## Died 18 0 0 6 124
## Euthanasia 236 0 25 200 704
## Return_to_owner 2088 0 5 880 617
## Transfer 1936 0 9 415 4695
## class.error
## Adoption 0.1152656
## Died 1.0000000
## Euthanasia 0.9785408
## Return_to_owner 0.7548747
## Transfer 0.3345145we see that even the random forest does not manage to predict Death and that, overall, there is still considerable class error, except for Adoption (which is nevertheless the most frequent occurence)
In terms of performances we have
pred = predict(model1rf, newdata = valid)
accrf = sum(pred == as.character(valid$OutcomeType))/length(pred)
accrf## [1] 0.6446974We see that the improvement over the decision tree is only marginal (at the level of the third decimal digit) indicating that perhaps all the information that could have been extracted from our aggregated variables has been used and that, due to the aggregation, the decision tree was not overfitting. Let’s check the variable importance
importance(model1rf)## MeanDecreaseGini
## AnimalType 376.16487
## AgeCat 656.35646
## IsPureBreed 36.75881
## operated 2038.86619
## season 117.46951
## daytype 99.36731
## time 398.11176As a last informative effort I want to see if, given the importance of the related aggreated variables, the Date/Time aggregation caused any information loss
formula2<- OutcomeType~ AnimalType + AgeCat + IsPureBreed + operated +
month + day + time
model2rf <- randomForest(formula2, data=train)which definitely seems to be the case!
pred2 = predict(model2rf, newdata = valid)
acc2rf = sum(pred2 == as.character(valid$OutcomeType))/length(pred2)
acc2rf## [1] 0.6509886importance(model2rf)## MeanDecreaseGini
## AnimalType 375.06582
## AgeCat 714.37886
## IsPureBreed 43.52733
## operated 2112.05688
## month 214.82721
## day 249.11009
## time 422.98671In order to try to understand a bit better how the day and month variables influence the outcome, let’s run a decision tree and visualize it
modeldt2 <-rpart(formula2, method="class", data=train)
fancyRpartPlot(modeldt2)
So, the story we learned so far is the following: if an animal has not been operated then most likely it gets transferred. If she has been operated then, if she is still a kitten/puppy, then adoption is most likely. If she’s older then the animal type enters into play. If she’s a dog then age matters as younger animals are more likely to end up adopted and older ones returned to their owners. If she’s a kitty then apparently what matters most is the time of the day. If we’re outside working hours then adoption is most likely (perhaps because people tend to visit shelters in their free time) whereas during working hours transfer is most likely.
No need to say, death and euthanasia go undetected, as well as return to owner if the animal is a cat. Also, the split on the hours suggests that a binary distinction between working hours vs non-working hours might be best.
Breed purity does not seem to be extremely relevant but this only tells us that purity is not the most significant breed feature and I am pretty sure that other extremely informative ‘macro’ features can be extracted and used.. So, more work on the breeds needs to be done!
# recast the final prediction in the format required for submission
model2rf <- randomForest(formula2, data=data)
prediction <- predict(model2rf, test, type = 'vote')
submission <- data.frame('ID' = test$ID, prediction)
write.csv(submission, 'rf2submission.csv', row.names = F)Any feedback is more than welcome!