startup_project.Rmd

---
title: "Startups Analysis"
author: "Iryna Popovych, Sofiya Hevorhyan"
date: 'April 2019'
output: html_notebook
---
<style type="text/css">

body{ /* Normal  */
  font-size: 16px;
  color: black;
  font-family: "Times New Roman", Times, serif;
  background-image: url(https://png.pngtree.com/thumb_back/fw800/back_pic/04/21/71/295830892804e4c.jpg);
  background-position: center center;
  background-attachment: fixed;
  background-repeat: no-repeat;
  background-size: 100% 100%;
  }
h1.title {
  font-size: 34px;
  color: DarkRed;
}
h1 { /* Header 1 */
  font-size: 28px;
  color: black;
}
h2 { /* Header 2 */
    font-size: 22px;
  font-family: "Times New Roman", Times, serif;
  color: darkpurple;
}
h3 { /* Header 3 */
  font-size: 18px;
  font-family: "Times New Roman", Times, serif;
  color: darkblue;
}
</style>


*This is the final project of Econometrics course taught by Liubomyr Bregman in Spring semester, 2019.*
*The work presented below was done by Sofiya Hevorhyan and Iryna Popovych, students of IT&BA program at UCU.*

# About
In this project we will work with the startup data. We have a dataset with information about 472 startups each having 116 characteristics: numerical, text, boolean etc. This data represents startups in terms of their industries, historical data about company, some details on history of investments, activities of companies, details about employees and technologies used.

Our aim is to practice our skills in data exploration, data cleaning and vizualization, feature creation and feature selection, treating missing values using the best practices, do some hypothesis testing, and, finally, build a model predict sucess of the company based on the information available.

Here we go!

## 1 Dealing with data
### The idea of the structure of the data 

```{r}
startup <- read.csv(file="./data/CAX_Startup_Data.csv", header=TRUE,as.is=TRUE)
head(startup)
```
### Function for changing character to two-stage factor
```{r}
# function for transforming character col -> two-stage factor col
my.col.change.as.factor <- function(dataframe, col.name) {
  # change all data in col to lowecase, without blank spaces
  dataframe[,col.name] <- trimws(tolower(dataframe[,col.name])) 
  vars <- unique(dataframe[,col.name])
  if (length(vars) == 2 || sum(is.na(vars)) == 1) {
    dataframe[,col.name] <- as.factor(dataframe[,col.name])
    return(dataframe)
  }
}
```
### Replacing common unknowns
We can see that some missing values are represented as "No Info" or are just "blanks". So, the first step is to replace "No Info" and "blanks" with "NA" to get better sense of missing values in data.
However, it doesn't mean that there are no others nontypical strings in other columns. We will check this later by unique() function
```{r}
startup[startup=="No Info"]<- NA 
startup[startup==""]<- NA 
startup[startup=="unknown amount"]<- NA 
```
### Retyping and renaming. Data convertation
```{r}
#Changing company_name (id) representation

library('stringr')

names(startup)[1] <- "id"
startup$id <- as.numeric(str_replace(startup$id, 
                                     pattern="Company",
                                     replacement = ""))
#names
names(startup) <- tolower(names(startup))
my.names <- names(startup)

#Changing dependent company status
startup <- my.col.change.as.factor(startup, my.names[2])
```
Now, instead of doing all these transformations by hand, let's just try to automate this process.
After analyzing variables in this data, we consider 4 types of features that belong to class 'character'
1) the string itself (e.g 'description'), that should remain so 
2) the string that should be converted to 2- or multi-stage factor (e.g. company status which is 2-stage)
3) the string that should be converted to numeric data type (e.g. year of founding)
4) the string that should be converted to date
```{r}
date.col <- c(13, 14)
date.col.names <- my.names[date.col]
startup$est..founding.date <- as.Date(startup$est..founding.date, "%m/%d/%Y") 
startup$last.funding.date <- as.Date(startup$last.funding.date, "%m/%d/%Y")

# it's hard to get some use of the date type in out model
# so we extract the year from our last.funding.date
startup$last.funding.date <- as.numeric(format(startup$last.funding.date, '%Y'))
```
### Retyping col with numbers to numeric data type
As we can see, there are no warnings, so, all the data is in rows are numbers
```{r}
num.col<- c(3:5, 10:11,15,
            18:23,25,61,66,68:70,72,74,88,92,94:96,98,99,102:116) 
num.col.names <- my.names[num.col]

for(i in num.col)
  {
   startup[,i]<-round(as.numeric(startup[,i]), 4) # round all numerical data to have clearer representation
}
```
### Retyping col to two-stage factor
Get all the columns that are supposed to have only 2 types of answers (Y/N questions mostly)
```{r}
two.fact.col <- c(2, 12, 24, 27, 29:32,
                  34, 36, 38, 40:42, 44:53,
                  55, 58, 63:64, 77:78, 81:86, 89:91)
for(i in two.fact.col)
  {
   startup <- my.col.change.as.factor(startup, my.names[i])
}
```
### Retyping col to multi-stage factor
```{r}
# in col 37 there is some inconstistency in data representation
# basically both "not applicable" and "no" mean the same, so we just replace it
startup[,37] <- str_replace(startup[,37],
                            pattern="not applicable",
                            replacement="no")
multi.fact.col <- c(26, 28, 33, 35, 37, 39, 43, 54, 56,
                    57, 59, 60, 65, 67, 71, 73, 75, 76,
                    79, 80, 87, 93, 97, 100, 101)

# some of the multi-stage factors should be strongly ordered
ordered.levels <- c(26, 28, 43, 56, 57, 59, 67, 71, 73, 75,76, 79, 80, 87, 93)

my.replace <- function(startup, col.name, patt, replac) {
  startup[,col.name] <- str_replace(startup[,col.name], 
                                      pattern=patt,
                                      replacement = replac)
  return(startup)
}

change.as.multi.factor <- function(startup, col.name) {
  "
  to lower different types of multi-stage factor (which we later should reorder) some of them are changed (small to low, etc.) to make reordering easier
  1) type: none, few, mane
  2) type: low, medium, high
  3) type: none, low, medium, high
  "
  startup[,col.name] <- trimws(tolower(startup[,col.name]))
  
  
  if (col.name %in% ordered.levels) {
    startup <- my.replace(startup, col.name, "nothing", "none")
    startup <- my.replace(startup, col.name, "small", "low")
    startup <- my.replace(startup, col.name, "large", "high")
    #little trick, by renaming only this phrase, the col now is right-ordered
    startup <- my.replace(startup, col.name, "average", "close_average")
  }
  
  
  vars <- unique(startup[,col.name])
  startup[,col.name] <- as.factor(startup[,col.name])
  col.levels <- levels(startup[, col.name])
  if (length(col.levels) == 3 && ("few" %in% col.levels)) {
    col.levels <- c("none", "few", "many")
  } else if (length(col.levels) == 3 && "low" %in% col.levels) {
    col.levels <- c("low", "medium", "high")
  } else if (length(col.levels) == 4 && "low" %in% col.levels) {
    col.levels <- c("none", "low", "medium", "high")
  }
  startup[,col.name] <- factor(startup[,col.name],
                               levels = col.levels)
  return(startup)
}

for(i in multi.fact.col)
  {
  startup <- change.as.multi.factor(startup, i)
}

#save text columns separately, now only these col are characters
char.col <- c(6:9, 16:17, 62) 
```
### Missing values
For modelling, we need relevant data with not much missing values. Let's take a look on how they are distributred here:
```{r}
#all incomplete rows
sum(!complete.cases(startup))
```
If we go by rows, most of the startups do not have ALL of the variables recorded. Let's determine percentage missing for each variable:
```{r}
# Percent missing value for each variable
mis_val<-sapply(startup, function(x) sum(is.na(x)))
percent_mis<-as.data.frame(round((mis_val/nrow(startup))*100,1))

name<-row.names(percent_mis)
pcnt_mis_var<-cbind(name,percent_mis)
row.names(pcnt_mis_var)<-NULL
colnames(pcnt_mis_var)<-c("variable","Percent.Missing")
```

### Separating data
Any variable missing more than 50% should not be used in modeling as it can give false impression of relationship with dependent and can pollute the model. We prefer to keep variable with less than 40% missing value only for treatment and anything above 40% missing values are to be used in testing only to give additional insights. 
```{r}
# keeping only variables with less than 40% missing
new_var<-as.character(pcnt_mis_var$variable[which(pcnt_mis_var$Percent.Missing<40)])
new_startup<-startup[new_var]
# separate data frame for more than 40% missing
other_var<-as.character(pcnt_mis_var$variable[which(pcnt_mis_var$Percent.Missing>=40)])
other_data<-startup[other_var]
```
We have separated 3 variables which are not suitable for further modelling.
```{r}
length(other_var)
```
Let's look how we reduced number of incomplete cases (rows which have missing values):
```{r}
sum(!complete.cases(startup)) - sum(!complete.cases(new_startup))
```
### Separate numerical and character data
It would be better if we separate out the numeric and character/factor variables from data frame. It would help in performing operations going forward. 
```{r}
# Separate data frame for numeric variables
cnt_var <- c(num.col.names, date.col.names)
cnt_var <- cnt_var[!cnt_var %in% other_var]
cnt_df<-new_startup[cnt_var]
# separate data frame for character variables
var <- colnames(new_startup) %in% cnt_var 
char_df <- new_startup[!var]
char_var <- names(char_df)
```
## Individual variables. Treating outliers
We need to remove ouliers or to subtract them with other values.
Take a look at quantiles for our numerical variables.
```{r}
library('matrixStats')
colquant <- as.data.frame(colQuantiles(data.matrix(cnt_df), probs = seq(0, 1, by= 0.1), na.rm = TRUE))
colquant
```
We create a function that replaces outliers with NA's using interquartile range rule. Outliers here are defined as observations that fall below Q1 - 1.5IQR or above Q3 + 1.5IQR.
```{r}
remove_outliers <- function(x, na.rm = TRUE, ...) {
  qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
  H <- 1.5 * IQR(x, na.rm = na.rm)
  y <- x
  y[x < (qnt[1] - H)] <- NA
  y[x > (qnt[2] + H)] <- NA
  return(y)
}
```
Now we loop through all the columns witth numeric entries and eliminate ouliers using our function:
```{r}
for (n in colnames(cnt_df)){
  if (class(cnt_df[,n][1]) != "Date"){
     col <- cnt_df[,n]
     summary(col)
     col <- remove_outliers(col)
     summary(col)
     cnt_df[,n] <- col
  }
}
```
Take a look on quantiles again, some of the upper quantiles, like 'Team.size.all.employees' have changed significatly, from 5000 to 110, or 'Last.Funding.Amount' from 7.7e+07 to 1.4e+07  - our job with ouliers is done!
```{r}
colquant <- as.data.frame(colQuantiles(data.matrix(cnt_df), probs = seq(0, 1, by= 0.1), na.rm = TRUE))
colquant
```

## Individual variables. Treating missing values
As we can see, every variable in our dataset has NA's.
```{r}
missing <- as.data.frame(sapply(cnt_df, function(x) sum(is.na(x))))
missing
```
Among of other methods, we have chosen to impute missing values using K-nearest-neighbours algorithm. k-Nearest Neighbour Imputation is based on a variation of the Gower Distance for numerical, categorical, ordered and semi-continous variables.
```{r}
library('VIM')
cnt_df <- kNN(cnt_df, imp_var = FALSE)
```
Fill missing dates by approximation, because VIM knn does not work with dates.
```{r}
#library('zoo')
#cnt[, 'Est..Founding.Date'] <- na.approx(cnt_df$Est..Founding.Date, na.rm = FALSE)
#cnt[, 'Last.Funding.Date '] <- na.approx(cnt_df$Last.Funding.Date, na.rm = FALSE)

#library(padr)
#pad(cnt_df, by='est..founding.date')

# --------------------------- TODO : DO SOMETHING WITH MISSING DATES ------------------------------

#print(cnt$est..founding.date)
```
We filled all the gaps! (EXCEPT DATES)
```{r}
missing <- as.data.frame(sapply(cnt_df, function(x) sum(is.na(x))))
missing
```
## Feature creation
From now on, we work with cnt_df table with clean and imputed data, with no outliers.
```{r}
# change scale to log
cnt_df$last.funding.amount <- log(cnt_df$last.funding.amount)

# Create additional features. For example, count number of investors for company
char_df$investor.count<-length(strsplit(char_df$investors, "|",fixed=T))
for (i in (1:length(char_df$investors))) {
  if(is.na(char_df$investors[i])==T){
    char_df$investor.count[i]<- NA}
  else{
    lst<-strsplit(char_df$investors[i], "|", fixed=T)
    char_df$investor.count[i]<-length(lst[[1]])
  }
}
# renew names of char_df
char_var <- names(char_df)

# lastfunding.amount vs. age of company ration
cnt_df$funding.agecomp.ratio <-
  cnt_df$last.funding.amount/cnt_df$age.of.company.in.years
# renew names of cnt_df
cnt_var <- names(cnt_df)
```
### PCA 
We decided to perform PCA to see what are the principal components. However, the dataset is too wide to provide good and valuable results which describe high percent of population variance. After dozens of failed attempts, we decided to perform PCA on small part of our data, in particular, just column characteristics that describe different skills percentage for founders and co-founders of startups.
```{r}
#head(cnt_df)

# cols 31-38 are percent skills 
cnt.pca <- prcomp(cnt_df[,c(27, 28, 29, 30, 31, 32, 33, 34, 35, 36)], center = TRUE,scale. = TRUE)

#cnt_df.pca <- prcomp(cnt[,c(3, 22, )], center = TRUE,scale. = TRUE)
summary(cnt.pca)
```
From the summary we can see that the first principal component is srill not very descriptive, because it describes only 28% of data variace. We can see that there are no clear groupd from the vizualization below, too:
```{r}
# install.packages('scales')
# library(devtools) 
# install_github("vqv/ggbiplot")

library(ggbiplot)

ggbiplot(cnt.pca)
```
PCA came out to be not very applicable in our case. After all, we tried to do ot and learned some new stuff while trying to do this.

## Graphics #1
```{r}
library(ggplot2)

# boxplot of employee count
boxplot(cnt_df$employee.count, main="box plot of employee count",
        ylab="Employee count")

# histogram with black outline, white fill and median line
ggplot(cnt_df, aes(x=employee.count))+
  geom_histogram(binwidth=5, colour="black", fill="white")+
  geom_vline(aes(xintercept=median(employee.count, na.rm=T)),
             color="red", linetype="dashed", size=1)+
  ggtitle("Histogram of Employee count")+
  xlab("Employee Count") +
  ylab("Frequency")+
  theme_light()
```
## Graphics #2
```{r}
# adding dependent variable to numeric data frame
cnt_df$dependent.company.status <-char_df$dependent.company.status
cnt_var <- names(cnt_df)

# box plot to see difference in mean of team size w.r.t two categories of dependent
ggplot(cnt_df, aes(x=dependent.company.status,y=team.size.all.employees,
                   fill=dependent.company.status)) +
  geom_boxplot()

# data preparation for bar chart
avg_emp<-aggregate(cnt_df$team.size.all.employees, 
                   by=list(cnt_df$dependent.company.status),
                   FUN=mean, na.rm=TRUE)
colnames(avg_emp)<-c("company.status","avg.employee.size")

# bar chart to check for difference in mean
ggplot(avg_emp, aes(x = company.status, y = avg.employee.size)) +
  geom_bar(stat = "identity")
```
## Graphics #3
```{r}
ggplot(cnt_df, aes(x=cnt_df$last.funding.amount)) + 
  geom_histogram(aes(fill=..count..), binwidth = 0.7) +
  scale_x_continuous(name="last funding amount, log",
                     breaks=seq(9, 17, by=1)) +
  scale_y_continuous(name="number of startups") +
  ggtitle("Distribution of last funding amount") +
  theme_bw() +
  geom_vline(xintercept = mean(cnt_df$last.funding.amount), size = 1,
             colour = "#FF3721",
             linetype = "dashed")

ggplot(cnt_df, aes(x = number.of.investors.in.seed,
                   fill=cnt_df$dependent.company.status)) +
  geom_histogram(aes(y=..count..), binwidth=1,
                 position="identity", alpha=0.5) +
  scale_x_continuous(name = "number of investors in seed") +
  scale_y_continuous(name = "number of startup") +
  ggtitle("Distribution of num of investors in seed with respect to success") +
  theme_bw() +
  scale_fill_brewer(palette="Pastel1")
```
## Graphics #4
```{r}
ggplot(cnt_df, aes(x=percent_skill_data.science)) +
  geom_histogram(fill="yellowgreen", color="yellow4",binwidth = 1.5) +
  facet_grid(dependent.company.status ~ .)+
  ggtitle("Data science skills influence on success")
```
## Variable Selection. Random Forest
Now, we are trying to take closer look at our numerical data to understand which of the variables are the most important. We use randomForest approach and caret varImp() function
```{r}
cnt_df$dependent.company.status <- as.numeric(cnt_df$dependent.company.status)-1

library(caret)
glm_model <- glm(formula = dependent.company.status ~ ., data = cnt_df)
importance1 <- varImp(glm_model)

library(randomForest)
rf_model = randomForest(dependent.company.status~., data=cnt_df[,-c(42,43)])
# Create an importance based on mean decreasing gini
importance2 <- importance(rf_model)

#ranks features, the bigger oberall, the more var is important
head(importance1)
head(importance2)
# open in large window
varImpPlot(rf_model)

# select features for modeling
name_var1 <- rownames(importance1)[apply(importance1, 1, function(x) x > 2)]
name_var2 <- rownames(importance2)[apply(importance2, 1, function(x) x > 3)]
mean(name_var1 != name_var2)
train1 <- cnt_df[, c(name_var1, "dependent.company.status")]
train2 <- cnt_df[, c(name_var2, "dependent.company.status")]

name_var <- c(name_var1, name_var2)
train.data <- cnt_df[, c(name_var, "dependent.company.status")]
```

## Statistical tests
```{r}
# t-test for checking difference in mean
t.test(team.size.all.employees~dependent.company.status, data=cnt_df)
# tabulating data for chi-sq test
tab<- table(char_df$dependent.company.status,char_df$local.or.global.player)
# chi-sq test for categorical variable
chisq.test(tab)
```
## Variable Importance

```{r}
library(devtools)
# install_github("tomasgreif/woe")
# install_github("tomasgreif/riv")
# install.packages("DBI",dependencies=TRUE)
# library(riv)
library(woe)
library(DBI)

cnt_df$fail <- as.numeric(!cnt_df$dependent.company.status)
# calculation of information value
IV<-iv.mult(cnt_df[, -c(42, 43, 45)],y="fail",summary=TRUE)

var<-IV[which(IV$InformationValue>0.1),]
var1<-var[which(var$InformationValue<0.5),]
var3<-var1$Variable
var3
# for model using information value
train3<-cnt_df[, c(var3, "dependent.company.status")]
```
## Model Building
```{r}
# fitting stepwise binary logistic regression with logit link function

mod1<-step(glm(dependent.company.status~., 
              family = binomial(link=logit),data = train.data))

mod2<-step(glm(dependent.company.status~., 
               family = binomial(link=logit),data = train3))

summary(mod1)
summary(mod2)

# final logistic regression model
model<-glm(formula = dependent.company.status ~ 
             number.of.investors.in.seed +
             team.size.senior.leadership + 
             number.of.investors.in.seed + 
             number.of.recognitions.for.founders.and.co.founders + 
             percent_skill_data.science + 
             percent_skill_business.strategy + 
             percent_skill_sales + 
             last.funding.date +
             internet.activity.score,
           family = binomial(link = logit), data = cnt_df)

# model summary
summary(model)
```
## Results
### Confusion matrix
```{r}
library(caret)

confusion.vs <- confusionMatrix(as.factor(round(mod1$fitted.values)),
                                     as.factor(cnt_df$dependent.company.status))


qplot(as.factor(cnt_df$dependent.company.status), 
      as.factor(round(mod1$fitted.values)), 
      colour= cnt_df$dependent.company.status, geom = c("boxplot", "jitter"), 
      main = "predicted vs. observed using VS", 
      xlab = "Observations", ylab = "Predictions") + 
  scale_color_gradientn(colors = c("red", "black"))

confusion.vi <- confusionMatrix(as.factor(round(mod2$fitted.values)),
                                     as.factor(cnt_df$dependent.company.status))

qplot(as.factor(cnt_df$dependent.company.status), 
      as.factor(round(mod2$fitted.values)), 
      colour= cnt_df$dependent.company.status, geom = c("boxplot", "jitter"), 
      main = "predicted vs. observed using VI", 
      xlab = "Observations", ylab = "Predictions") + 
  scale_color_gradientn(colors = c("purple", "black"))


confusion.vs
confusion.vi

```
```{r}
# model fit (Hosmer and Lemeshow goodness of fit (GOF) test)
library(ResourceSelection)
hoslem.test(cnt_df$dependent.company.status, mod1$fitted.values, g=10)
hoslem.test(cnt_df$dependent.company.status, mod2$fitted.values, g=10)
```