install.packages(c("psych", "polycor", "lmtest", "car"))
library(psych)
library(polycor)
library(lmtest)
library(car)
#-------------------------------------------------------------------------------
# DATASET LOADING:

# factorial analysis 
instrument <- read.csv2("G:\\Mi unidad\\1 trabajo\\3 Investigación\\2 Proyectos de Investigación finalizados\\Sin Fondos\\2025 Indice de Oferta Servicios Digitales Municipios\\MDOI\\PeerJ\\BASES DE DATOS y codigos\\Instrument.csv",dec=",")

# regression
municipalities <- read.csv2("G:\\Mi unidad\\1 trabajo\\3 Investigación\\2 Proyectos de Investigación finalizados\\Sin Fondos\\2025 Indice de Oferta Servicios Digitales Municipios\\MDOI\\PeerJ\\BASES DE DATOS y codigos\\Municipalities.csv",dec=",")

#-------------------------------------------------------------------------------
# DATASET PREPARATION

dim(instrument)
instrument <-instrument[,2:89]
summary(instrument)

municipalities$Region <- factor(municipalities$Region)
municipalities$Region <- relevel(municipalities$Region,ref="Metropolitana")
summary(municipalities)

#-------------------------------------------------------------------------------
# GRAPHS

# Scree Plot
scree(instrument, main = "Scree plot")

# Parallel Analysis
fa.parallel(instrument, fa = "fa", n.iter = 100, show.legend = TRUE, main = "Paralel analysis for EFA")


# Histogram
hist(municipalities$MDOI_Index, prob=T,main="",xlab="MDOI",ylab="Frequency")
lines(density(municipalities$MDOI_Index),lwd=2)

#Boxplot
par(mar=c(8, 4, 2, 1))
boxplot(municipalities$MDOI_Index~ municipalities$Region,xlab="", ylab="MDOI index",las=3,municipalities=municipalities)
boxplot(municipalities$MDOI_Index)

#-------------------------------------------------------------------------------
# FACTORIAL ANALYSIS (section 4.1)

# Step 2: Estimate the tetrachoric matrix
matriz_tetra <- tetrachoric(instrument)$rho


# Stage 3:  (EFA) con oblimin rotation
n_factores <- 15 # number of factors

efa <- fa(r = instrument,
          nfactors = n_factores,
          fm = "minres",         
          rotate = "oblimin",    
          cor = "tet",
          scores = "regression")

print(efa)

#-------------------------------------------------------------------------------
# MIFL AND MDOI INDEXES (section 4.2)

# only values greater than 4 are kept in the loadings matrix (FL)
FL <- as.matrix(unclass(efa$loadings))
for(i in 1:dim(FL)[1]){
  for(j in 1:dim(FL)[2]){
    if(FL[i,j]<0.4){FL[i,j]=0}
  }
}
# keep only the highest loading per variable in the loadings matrix (FL 
for(i in 1:dim(FL)[1]){
  for(j in 1:dim(FL)[2]){
    if(FL[i,j]<max(FL[i,])){FL[i,j]=0}
  }
}

# divide by the sum of loadings per factor (SFL)
idx.factor <- matrix(c(0),88,15)
for(i in 1:88){
  for(j in 1:15){
    if(FL[i,j]>0){idx.factor[i,j]=j}
  }
}
idx.factor <- apply(idx.factor,1,sum)
SFL <- rep(0,88)
for(j in 1:15){
  SFL[which(idx.factor==j)] <- sum(FL[which(idx.factor==j),j]) 
}

for(i in 1:88){
  for(j in 1:15){
    FL[i,j] <- FL[i,j]/SFL[j]
    
  }
}

# calculation of MIFL index's
MIFL    <- as.matrix(instrument) %*% FL
mean    <-apply(MIFL,2,mean)
StdDev  <- apply(MIFL,2,sd )
Max     <- apply(MIFL,2,max)
Min     <-apply(MIFL,2,min)
perc.min.03 <-rep(0,15)
for(j in 1:15){
  aux <- MIFL[,j]
  perc.min.03[j] <- length(aux[aux<0.25])/344*100
}

#summary of MIFL
Tabla3 <- round(data.frame(mean,StdDev,Max,Min,perc.min.03),2)
Tabla3
# calculation of MDOI index's
EV <- efa$Vaccounted[2,]
MDOI <- rep(0,344)
for(j in 1:15){
  MDOI <- MDOI + MIFL[,j]*EV[j]
}


#-------------------------------------------------------------------------------
# LINEAR REGRESION (section 4.3)

model <- lm(log(MDOI_Index) ~ log(CDI) + log(Size) + log(No_Schooling) + log(Income) + Region , data=municipalities)
summary(model)
bptest(model)

hccm(model) # Estimador de la varianza robusto de White
tabla <- cbind( coef(model),sqrt(diag(hccm(model))),coef(model)/sqrt(diag(hccm(model))) )
colnames(tabla) <- c("coefficient", "standard error","T")
tabla
# ANOVA
anova <- data.frame(as.vector(summary(aov(model))[[1]][,2]/sum(summary(aov(model))[[1]][,2])*100))
colnames(anova) <- "anova"
rownames(anova) <- c("log(CDI)", "log(Size)", "log(No_Schooling)", "log(Income) ", "Region", "Residuals")
anova
#-------------------------------------------------------------------------------