## read the data from the Kaggle example (Polish company data, 
## see Readme file)

train <- read.csv("bankruptcy_Train.csv")

## cleanup: eliminate 10 largest absolute values
#train <- train[-6424,]
for(i in 1:10) {
  train <- train[-unique(apply(train,2,
                               FUN=function(x) which.max(abs(x)))
  ),]
}

## Recupero le variabili "di Altman" dal dataset

## Indice di flessibilità (Attivo corrente/Totale attivo)
train$x1 <- train$Attr3

## Indice di autofinanziamento (Utile non distribuito/totale attivo)
train$x2 <- train$Attr6

## ROA
train$x3 <- train$Attr7

## Patrimonio netto/Passività totali
train$x4 <- train$Attr8

## Vendite nette/Totale attivo
train$x5 <- train$Attr9

## Equazione di Altman
train$z.altman <- 1.2*train$x1 + 1.4*train$x2 + 3.3*train$x3 +
  0.6*train$x4 + 0.999*train$x5

summary(train$z.altman)
sum(train$z.altman < 1.81)/length(train$z.altman)*100

## Equazione di Altman, mercati emergenti
train$z.altman.2 <- 3.25 + 6.56*train$x1 + 3.26*train$x2 +
  6.72*train$x3 + 1.05*train$x4

summary(train$z.altman.2)
sum(train$z.altman.2 < 1.1)/length(train$z.altman.2)*100

alt.data <- train[, c("class", paste("x", 1:5, sep=""))]

## Variable summary (cfr. Altman 2000, Table 1)
## NB: x1-4 non sono espresse in percentuale per omogeneità
## con la formulazione di Biggeri et al., Es. 8.8
round(apply(alt.data, 2,
            function(x) tapply(x, alt.data$class, mean)), 3)

## Estrazione dei due gruppi:
## numerosità delle insolventi
n.1 <- sum(train$class)
insolventi <- train[train$class==1,]
sane.all <- train[train$class==0,]
## campione casuale di sane, stesso numero
sane <- sane.all[sample(1:dim(sane.all)[[1]], n.1), ]
adset <- rbind(insolventi, sane)

## AD lineare (Linear Discriminant Analysis, o LDA)
library(MASS)
adl <- lda(class ~ x1 + x2 + x3 + x4 + x5, data=adset,
           prior=c(1,1)/2)
adl$prior
adl$counts
adl$means

adl$scaling

## previsione del risultato:
padl0 <- predict(adl)

## risultato della classificazione:
head(padl0$class, 10)

## probabilità a posteriori:
head(padl0$posterior, 10)


## Gruppi di dimensioni diseguali: teniamo conto della
## diversa probabilità a priori di essere 0 o 1 mediante
## il Th. di Bayes (v. 8.6.4, Eq. 25)
adl.b <- lda(class ~ x1 + x2 + x3 + x4 + x5, data=train,
             prior=c(sum(train$class==0)/length(train$class),
                     sum(train$class==1)/length(train$class)))


####### PCA #######

## PCA example, kaggle data

## Advanced PCA with FactoMineR

library(FactoMineR)
pl.pca1 <- PCA(train[, paste("Attr",1:64, sep="")])

## plot map of individuals
plot(pl.pca1, select=1)

## inspect results
library(factoextra)
get_pca_ind(pl.pca1)$coord[1:10,]
get_pca_var(pl.pca1)$coord[1:10,]



## analysis of individuals' positioning  
fviz_pca_ind(pl.pca1, geom.ind = "point", pointshape = 21, 
             pointsize = 2, 
             fill.ind = cut(train$class,2), 
             col.ind = "black", 
             palette = "jco", 
             addEllipses = TRUE,
             label = "var",
             col.var = "black",
             repel = TRUE,
             legend.title = "Diagnosis") +
  ggtitle("2D PCA-plot from 30 feature dataset") +
  theme(plot.title = element_text(hjust = 0.5))

## analysis of variables
fviz_pca_var(pl.pca1, geom.ind = "point", pointshape = 21, 
             pointsize = 2, 
             fill.ind = train$class, 
             col.ind = "black", 
             palette = "jco", 
             addEllipses = TRUE,
             label = "var",
             col.var = "black",
             repel = TRUE,
             legend.title = "Diagnosis")

###### linear probability model #######

fm <- class ~ x1 + x2 + x3 + x4 + x5

lpm <- lm(fm, data=train)
summary(lpm)

head(fitted(lpm))

## range di y.hat
summary(fitted(lpm))

####### logit ########

logitm <- glm(fm, data=train, family=binomial(link="logit"))
summary(logitm)

## range di y.hat
summary(fitted(logitm))

## McFadden's pseudo-R2
logit0 <- glm(update(fm, .~1), data=train,
              family=binomial(link="logit"))

pR2 <- 1-logLik(logitm)/logLik(logit0)

##### probit ######

probitm <- glm(fm, data=train, family=binomial(link="probit"))
summary(probitm)

## confronto modelli:
library(texreg)
screenreg(list("logit"=logitm, "probit"=probitm, "LPM"=lpm))



################ glm (logit) on full data ############

## extract first k prin comps
pl.pc <- pl.pca1$ind$coord

logitmod <- glm(train$class ~ pl.pc, family=binomial(link="logit"))
summary(logitmod)

## logit with 64 variables:
logitmod64 <- glm(train$class ~ as.matrix(train[,-65]), family=binomial(link="logit"))
summary(logitmod64)