# R codes used in the thesis
# Logistic regression for the whole data set

setwd("D:/Thesis/Data in Tunisia")
Tunisia.data <- read.csv("data.csv", head=TRUE)
attach(Tunisia.data)
save( Tunisia.data, file="Tunisia.data.rda" )
data.logit <- glm(Y ~ AGE + GENDER + EDUC + CIVSTATUS + DEPCHILD + TYPBORR
                  + TYPCONTR + OBJLOAN + SOCILEVEL + IMPROVEMENT + PROBLEM
                  + REPAYMENT + KINDIMF + USEMICRO + FINAINCLUS + SAVING
                  + USESAVING + COLLATERAL + OTHERLOANS + INDGROUP + BUSISECTOR
                  + REA.ACTIVITY + REA.ASKLOAN, data = Tunisia.data, family = binomial())
summary(data.logit)

# Commpute confident interval of coefficients, OR and confidence interval of OR in the whole data
round(cbind(confint.default(data.logit)),3)
round(cbind(exp(cbind(OR=data.logit$coefficients)), exp(confint.default(data.logit))),3)

# AIC Backward Stepwise Elimination
# Running "stepAIC" for the whole data

library(MASS)
dataAIC.step <- stepAIC(data.logit, trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary(dataAIC.step)
AICmodel <- glm(Y ~ GENDER + EDUC + TYPCONTR +OBJLOAN + PROBLEM
                + REPAYMENT + KINDIMF + FINAINCLUS + SAVING +USESAVING
                + COLLATERAL +INDGROUP + BUSISECTOR, data = Tunisia.data,
                family = binomial())
round(cbind(confint.default(dataAIC.step)),3)
round(cbind(exp(cbind(OR=dataAIC.step$coefficients)), exp(confint.default(dataAIC.step))),3)
# Delete variable in AIC optimal model by Backward stepwise elimination procedure

# Delete "USESAVING", p-value = 0.182
USESAVING.logit <- glm(Y ~ GENDER + EDUC + TYPCONTR + OBJLOAN + PROBLEM
                       + REPAYMENT + KINDIMF + FINAINCLUS + SAVING + COLLATERAL
                       + INDGROUP + BUSISECTOR, data = Tunisia.data, family = binomial())
USESAVING.AIC <- stepAIC(USESAVING.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary(USESAVING.AIC)

# Delete "REPAYMENT", p-value = 0.148
REPAYMENT.logit <- glm(Y ~ GENDER + EDUC + TYPCONTR +OBJLOAN + PROBLEM
                       + KINDIMF + FINAINCLUS + SAVING + COLLATERAL +INDGROUP
                       + BUSISECTOR, data = Tunisia.data, family = binomial())
REPAYMENT.AIC <- stepAIC(REPAYMENT.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary(REPAYMENT.AIC)

# Delete "PROBLEM", p-value = 0.262
PROBLEM.logit <- glm(Y ~ GENDER + EDUC + TYPCONTR +OBJLOAN + KINDIMF
                     + FINAINCLUS + SAVING + COLLATERAL +INDGROUP + BUSISECTOR,
                     data = Tunisia.data, family = binomial())
PROBLEM.AIC <- stepAIC(PROBLEM.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary(PROBLEM.AIC)

# Delete "BUSISECTOR", p-value = 0.095
BUSISECTOR.logit <- glm(Y ~ GENDER + EDUC + TYPCONTR +OBJLOAN + KINDIMF
                        + FINAINCLUS + SAVING + COLLATERAL +INDGROUP,
                        data = Tunisia.data, family = binomial())
BUSISECTOR.AIC <- stepAIC(BUSISECTOR.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary(BUSISECTOR.AIC)

# Delete "INDGROUP ", p-value = 0.124
INDGROUP.logit <- glm(Y ~ GENDER + EDUC + TYPCONTR +OBJLOAN + KINDIMF
                      + FINAINCLUS + SAVING + COLLATERAL, data = Tunisia.data, family = binomial())
INDGROUP.AIC <- stepAIC(INDGROUP.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary(INDGROUP.AIC)

# Delete "TYPCONTR", p-value = 0.11
TYPCONTR.logit <- glm(Y ~ GENDER + EDUC +OBJLOAN + KINDIMF + FINAINCLUS 
                      + SAVING + COLLATERAL , data = Tunisia.data, family = binomial())
TYPCONTR.AIC <- stepAIC( TYPCONTR.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary( TYPCONTR.AIC)

# Delete "COLLATERAL", p-value = 0.119
COLLATERAL.logit <- glm(Y ~ GENDER + EDUC +OBJLOAN + KINDIMF + FINAINCLUS +
                          SAVING, data = Tunisia.data, family = binomial())
COLLATERAL.AIC <- stepAIC( COLLATERAL.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary( COLLATERAL.AIC)

# Delete "KINDIMF", p-value = 0.051
KINDIMF.logit <- glm(Y ~ GENDER + EDUC +OBJLOAN + FINAINCLUS + SAVING, data = Tunisia.data, family = binomial())
KINDIMF.AIC <- stepAIC(KINDIMF.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary(KINDIMF.AIC)

# Delete "OBJLOAN ", p-value = 0.014
OBJLOAN.logit <- glm(Y ~ GENDER + EDUC + FINAINCLUS + SAVING , data = Tunisia.data, family = binomial())
OBJLOAN.AIC <- stepAIC(OBJLOAN.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary(OBJLOAN.AIC)

# Delete "GENDER", p-value = 0.023
GENDER.logit <- glm(Y ~ EDUC + FINAINCLUS + SAVING , data = Tunisia.data, family = binomial())
GENDER.AIC <- stepAIC(GENDER.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary( GENDER.AIC)

#Delete "EDUC ", p-value = 0.0048
EDUC.logit <- glm(Y ~ FINAINCLUS + SAVING, data = Tunisia.data, family = binomial())
EDUC.AIC <- stepAIC(EDUC.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary( EDUC.AIC)

# Delete "FINAINCLUS ", p-value = 0.979
FINAINCLUS.logit <- glm(Y ~ SAVING, data = Tunisia.data, family = binomial())
FINAINCLUS.AIC <- stepAIC(FINAINCLUS.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary( FINAINCLUS.AIC)

# Delete "SAVING", p-value < 2e-16
SAVING.logit <- glm(Y ~ 1 , data = Tunisia.data, family = binomial())
SAVING.AIC <- stepAIC(SAVING.logit , trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary(SAVING.AIC)

# BIC Backward Stepwise Elimination
# Running "stepAIC" for the whole data, "k = log (n)"

library(MASS)
dataBIC.step <- stepAIC(data.logit, trace = TRUE, k=log(length(Y)), data=Tunisia.data, direction="both")
summary(dataBIC.step)
round(cbind(confint.default(dataBIC.step)), 3 ); round(cbind(exp(cbind(OR=dataBIC.step$coefficients)),
                                                             exp(confint.default(dataBIC.step))),3)

# Delete variable in BIC optimal model by Backward stepwise elimination procedure

#Delete "OBJLOAN ", p-value = 0.0138
OBJLOAN.logit <- glm(Y ~ GENDER + EDUC + FINAINCLUS + SAVING , data = Tunisia.data, family = binomial())
OBJLOAN.BIC <- stepAIC(OBJLOAN.logit , trace = TRUE, k=log(length(Y)), data=Tunisia.data, direction="both")
summary(OBJLOAN.BIC)

# Delete "GENDER", p-value = 0.0226
GENDER.logit <- glm(Y ~ EDUC + FINAINCLUS + SAVING , data = Tunisia.data, family = binomial())
GENDER.BIC <- stepAIC(GENDER.logit, trace = TRUE, k=log(length(Y)), data=Tunisia.data, direction="both")
summary( GENDER.BIC)

# Delete "EDUC", p-value = 0.0048
EDUC.logit <- glm(Y ~ FINAINCLUS + SAVING, data = Tunisia.data, family = binomial())
EDUC.BIC <- stepAIC(EDUC.logit , trace = TRUE, k=log(length(Y)), data=Tunisia.data, direction="both")
summary( EDUC.BIC)

# Delete "FINAINCLUS", p-value = 0.979
FINAINCLUS.logit <- glm(Y ~ SAVING, data = Tunisia.data, family = binomial())
FINAINCLUS.BIC <- stepAIC(FINAINCLUS.logit , trace = TRUE, k= log(length(Y)), data=Tunisia.data, direction="both")
summary( FINAINCLUS.BIC)

# Delete "SAVING ", p-value < 2e-16
SAVING.logit <- glm(Y ~ 1 , data = Tunisia.data, family = binomial())
SAVING.BIC <- stepAIC(SAVING.logit , trace = TRUE, k= log(length(Y)), data=Tunisia.data, direction="both")
summary(SAVING.BIC)

# Plot AIC values of models obtained from deleting or adding variables to the full model
numberofvar <- c(23:0)
AICvalue <- c(309.12,307.12,305.15,303.22,301.41,299.62,297.82 ,296.79,294.98,294.27,293.88,
              293.89, 293.90,293.45,294.27,294.62,295.14,295.66,297.34,301.49,305.06 ,311.51,336.40,519.48)
plot(numberofvar, AICvalue, col ="blue",pch = 20, xlab="Number of variables in the models",
     ylab="AIC of coresponding models")
lines(numberofvar, AICvalue)

# Plot BIC values of models obtained from deleting or adding variables to the full model
numberofvar <- c(23:0)
BICvalue <- c(390.46,385.07,379.71, 374.39,369.19,364.02,358.82,354.41,349.20,345.11, 341.33,337.95,
              334.56,330.73,328.16, 325.12,322.25,319.38, 317.68,318.43,318.62,321.68,343.17,522.87)
plot(numberofvar,BIC, col ="blue",pch = 20, xlab="Number of variables in the models", ylab="AIC
     of coresponding models")
lines(numberofvar,BICvalue)

# Create 30 sub-samples to compute P{Y=1|X} to test the fitness of the model and 
# test the stability of variables in AIC optimal model
u=1:404;
Se.learning=sample(u,300, replace=FALSE)
learning.data=Tunisia.data[Se.learning,]
remain=u[-Se.learning]
test.data=Tunisia.data[remain,]
learning.logit <- glm( Y ~ AGE + GENDER + EDUC + CIVSTATUS + DEPCHILD + TYPBORR
                       + TYPCONTR + OBJLOAN + SOCILEVEL + IMPROVEMENT + PROBLEM
                       + REPAYMENT + KINDIMF + USEMICRO + FINAINCLUS + SAVING
                       + USESAVING + COLLATERAL + OTHERLOANS + INDGROUP + BUSISECTOR
                       + REA.ACTIVITY + REA.ASKLOAN, data = learning.data, family = binomial())
learning.stepAIC= stepAIC(learning.logit, trace = 1, keep = NULL, k=2, data=learning.data, direction="both");
summary(learning.stepAIC);
Allvar=names(data.logit$coefficients)
Selectvar=names(learning.stepAIC$coefficients);
SelectBeta=learning.stepAIC$coefficients;
Res<-numeric(length(Allvar))
for (i in 1 : length(Allvar) )
  for (j in 1: length(Selectvar))
    if (Allvar[i] == Selectvar[j])
      Res[i]<- SelectBeta[j]
Res[i]<- 0
Beta=as.matrix(Res)
const=matrix(data=1, nrow=length(remain), ncol=1)

test <- cbind(const,test.data$AGE, test.data$GENDER,test.data$EDUC,test.data$CIVSTATUS,
              test.data$DEPCHILD,test.data$TYPBORR, test.data$TYPCONTR,test.data$OBJLOAN,
              test.data$SOCILEVEL, test.data$IMPROVEMENT, test.data$PROBLEM,test.data$REPAYMENT,
              test.data$KINDIMF,test.data$USEMICRO,test.data$FINAINCLUS , test.data$SAVING,
              test.data$USESAVING,test.data$COLLATERAL, test.data$OTHERLOANS, test.data$INDGROUP,
              test.data$BUSISECTOR,test.data$REA.ACTIVITY, test.data$REA.ASKLOAN)
Xt=as.matrix(test)   # (input matrix of test sample)
Yt=as.matrix(test.data$Y) # (output vector of test sample)

# Using Beta obtained from learning.step to compute p(x) = P{Y = 1|X} and the number of OK pairs
pi.test=exp(Xt%???%Beta)/(1+exp(Xt%*%Beta))
a <- 0; b <- 0
for (i in 1 : length(row(Yt)))
  if (Yt[i] == 1 & pi.test[i] >= 0.5)
    a = a+1
for (i in 1 : length(row(Yt)))
  if (Yt[i] == 0 & pi.test[i] < 0.5)
    b= b+1
c <- 0; d <- 0
for (i in 1 : length(row(Yt)))
  if (Yt[i] == 1 & pi.test[i] < 0.5)
    c = c+1
for (i in 1 : length(row(Yt)))
  if (Yt[i] == 0 & pi.test[i] >= 0.5)
    d= d+1
a+b
c+d

# Compute pi for the model obtained fromthe whole data using AIC criterion
dataAIC= stepAIC(data.logit, trace = 1, keep = NULL, k=2, data=Tunisia.data, direction="both")
summary(dataAIC)
Allvar=names(data.logit$coefficients)
Selectvar=names(dataAIC$coefficients);
SelectBeta=dataAIC$coefficients;
Res <- numeric(length(Allvar))
for (i in 1 : length(Allvar) )
  for (j in 1: length(Selectvar))
    if (Allvar[i] == Selectvar[j])
      Res[i]<- SelectBeta[j]
Res[i] <- 0
Beta=as.matrix(Res)
const=matrix(data=1, nrow=404, ncol=1)
whole <- cbind(const,Tunisia.data$AGE, Tunisia.data$GENDER,Tunisia.data$EDUC,Tunisia.data$CIVSTATUS, 
               Tunisia.data$DEPCHILD,Tunisia.data$TYPBORR,Tunisia.data$TYPCONTR, Tunisia.data$OBJLOAN, 
               Tunisia.data$SOCILEVEL,Tunisia.data$IMPROVEMENT, Tunisia.data$PROBLEM, Tunisia.data$REPAYMENT,
               Tunisia.data$KINDIMF,Tunisia.data$USEMICRO,Tunisia.data$FINAINCLUS,Tunisia.data$SAVING, 
               Tunisia.data$USESAVING,Tunisia.data$COLLATERAL,Tunisia.data$OTHERLOANS, Tunisia.data$INDGROUP, 
               Tunisia.data$BUSISECTOR,Tunisia.data$REA.ACTIVITY, Tunisia.data$REA.ASKLOAN)
X=as.matrix(whole) # (input matrix of whole data)
Y=as.matrix(Tunisia.data$Y) # (output vector of whole data)
pi=exp(X%???%Beta)/(1+exp(X%*%Beta))
a <- 0; b <- 0
for (i in 1 : length(row(Y)))
  if (Y[i] == 1 & pi[i] >= 0.5)
    a = a+1
for (i in 1 : length(row(Y)))
  if (Y[i] == 0 & pi[i] < 0.5)
    b= b+1
c <- 0; d <- 0
for (i in 1 : length(row(Y)))
  if (Y[i] == 1 & pi[i] < 0.5)
    c = c+1
for (i in 1 : length(row(Y)))
  if (Y[i] == 0 & pi[i] >= 0.5)
    d= d+1
a+b
c+d

# Create 30 sub-samples to compute P{Y=1|X} to test the fitness of the model and test
#   the stability of variables in BIC optimal model
u=1:404
Se.learning=sample(u,300, replace = FALSE)
learning.data=Tunisia.data[Se.learning,]
remain=u[-Se.learning]
test.data=Tunisia.data[remain,]

learning.logit <- glm( Y ~ AGE + GENDER + EDUC + CIVSTATUS + DEPCHILD + TYPBORR
                       + TYPCONTR + OBJLOAN + SOCILEVEL + IMPROVEMENT + PROBLEM + REPAYMENT
                       + KINDIMF + USEMICRO + FINAINCLUS + SAVING + USESAVING + COLLATERAL
                       + OTHERLOANS + INDGROUP + BUSISECTOR + REA.ACTIVITY + REA.ASKLOAN,
                       data = learning.data, family = binomial())
learning.stepBIC= stepAIC(learning.logit, trace = TRUE, k=log(300), data=learning.data, direction="both")
summary(learning.stepBIC)
Allvar=names(data.logit$coefficients)
BICSelectvar=names(learning.stepBIC$coefficients);
BICSelectBeta=learning.stepBIC$coefficients;
Res <- numeric(length(Allvar))
for (i in 1 : length(Allvar) )
  for (j in 1: length(BICSelectvar))
    if (Allvar[i] == BICSelectvar[j])
      Res[i] <- BICSelectBeta[j]
Res[i] <- 0
BICbeta=as.matrix(Res)
const=matrix(data=1, nrow=length(remain), ncol=1)

test <- cbind(const,test.data$AGE, test.data$GENDER,test.data$EDUC,test.data$CIVSTATUS,
              test.data$DEPCHILD,test.data$TYPBORR, test.data$TYPCONTR,test.data$OBJLOAN,
              test.data$SOCILEVEL, test.data$IMPROVEMENT, test.data$PROBLEM,test.data$REPAYMENT,
              test.data$KINDIMF,test.data$USEMICRO,test.data$FINAINCLUS , test.data$SAVING,
              test.data$USESAVING,test.data$COLLATERAL, test.data$OTHERLOANS, test.data$INDGROUP,
              test.data$BUSISECTOR,test.data$REA.ACTIVITY, test.data$REA.ASKLOAN)
Xt=as.matrix(test) #(input matrix of test sample)
Yt=as.matrix(test.data$Y) # (output vector of test sample)

# Using Beta obtained from learning.stepBIC to compute p(x) = P{Y = 1|X} and the number of OK pairs

BICpi.test=exp(Xt%???%BICbeta)/(1+exp(Xt%???%BICbeta))
a <- 0; b <- 0
for (i in 1 : length(row(Yt)))
  if (Yt[i] == 1 & BICpi.test[i] >= 0.5)
    a = a+1
for (i in 1 : length(row(Yt)))
  if (Yt[i] == 0 & BICpi.test[i] < 0.5)
    b= b+1
c <- 0; d <- 0
for (i in 1 : length(row(Yt)))
  if (Yt[i] == 1 & BICpi.test[i] < 0.5)
    c = c+1
for (i in 1 : length(row(Yt)))
  if (Yt[i] == 0 & BICpi.test[i] >= 0.5)
    d= d+1
a+b
c+d

# Compute pi for the model obtained from the whole data using BIC criterion

dataBIC= stepAIC(data.logit, trace = TRUE, k=log(length(Y)), data=Tunisia.data, direction="both");
summary(dataBIC)
Allvar=names(data.logit$coefficients)
Selectvar=names(dataBIC$coefficients);
SelectBeta=dataBIC$coefficients;

Res <- numeric(length(Allvar))
for (i in 1 : length(Allvar) )
  for (j in 1: length(Selectvar))
    if (Allvar[i] == Selectvar[j])
      Res[i] <- SelectBeta[j]
Res[i] <- 0
Beta=as.matrix(Res)

const=matrix(data=1, nrow=404, ncol=1)
whole <- cbind(const,Tunisia.data$AGE, Tunisia.data$GENDER,Tunisia.data$EDUC,Tunisia.data$CIVSTATUS, 
               Tunisia.data$DEPCHILD,Tunisia.data$TYPBORR, Tunisia.data$TYPCONTR, Tunisia.data$OBJLOAN, 
               Tunisia.data$SOCILEVEL, Tunisia.data$IMPROVEMENT, Tunisia.data$PROBLEM, Tunisia.data$REPAYMENT,
               Tunisia.data$KINDIMF,Tunisia.data$USEMICRO,Tunisia.data$FINAINCLUS,Tunisia.data$SAVING, 
               Tunisia.data$USESAVING,Tunisia.data$COLLATERAL,Tunisia.data$OTHERLOANS, Tunisia.data$INDGROUP, 
               Tunisia.data$BUSISECTOR, Tunisia.data$REA.ACTIVITY, Tunisia.data$REA.ASKLOAN) 
X=as.matrix(whole) # (input matrix of whole data)
Y=as.matrix(Tunisia.data$Y) # (output vector of whole data)
pi=exp(X%???%Beta)/(1+exp(X%???%Beta))
a <- 0; b <- 0
for (i in 1 : length(row(Y)))
  if (Y[i] == 1 & pi[i] >= 0.5)
    a = a+1
for (i in 1 : length(row(Y)))
  if (Y[i] == 0 & pi[i] < 0.5)
    b= b+1
c <- 0; d <- 0
for (i in 1 : length(row(Y)))
  if (Y[i] == 1 & pi[i] < 0.5)
    c = c+1
for (i in 1 : length(row(Y)))
  if (Y[i] == 0 & pi[i] >= 0.5)
    d= d+1
a+b
c+d

# Plot the frequency of variables in 30 AIC samples
AICsample <- matrix(c(2,26,30,4,2,4,12,28,3,3,12,17,19,2,30,29,13,23,0,15,14,5,12),nrow = 1,
                    ncol = 23,byrow=TRUE)
barplot(AICsample, main = "Frequency of variables appearing in 30 AIC learning optimal models",
        beside = TRUE, names.arg =c("AGE", "GENDER", "EDUC", "CIVSTATUS ",
                                    "DEPCHILD", "TYPBORR", "TYPCONTR","OBJLOAN", "SOCILEVEL", "IMPROVEMENT",
                                    "PROBLEM", "REPAYMENT", "KINDIMF", "USEMICRO", "FINAINCLUS ",
                                    "SAVING", "USESAVING", "COLLATERAL", "OTHERLOANS", "INDGROUP",
                                    "BUSISECTOR", "REA.ACTIVITY ", "REA.ASKLOAN"), xlab="Variables", ylab="Number of variables")

# Plot the frequency of variables in 30 BIC samples
BICsample <- matrix(c(0,13,27,0,0,0,0,10,0,0,0,2,3,0,30,30,6,2,0,0,0,0,7),nrow = 1,
                    ncol = 23,byrow=TRUE)
barplot(BICsample, main = "Frequency of variables appearing in 30 BIC learning optimal models",
        beside = TRUE, names.arg =c("AGE", "GENDER", "EDUC", "CIVSTATUS ",
                                    "DEPCHILD", "TYPBORR", "TYPCONTR","OBJLOAN", "SOCILEVEL", "IMPROVEMENT",
                                    "PROBLEM", "REPAYMENT", "KINDIMF", "USEMICRO", "FINAINCLUS ",
                                    "SAVING", "USESAVING", "COLLATERAL", "OTHERLOANS", "INDGROUP",
                                    "BUSISECTOR", "REA.ACTIVITY ", "REA.ASKLOAN"), xlab="Variables", ylab="Number of variables")