Repaso

Resolvamos las siguientes preguntas:

1. Si un investigador usa la data de The Economist e hipotetiza que la calidad del proceso electoral depende del nivel de funcionamiento de gobierno y el nivel de participación política. ¿Ud qué puede afirmar?

rm(list = ls()) # clean memory
LINK='https://en.wikipedia.org/wiki/Democracy_Index'
PATH='//*[@id="mw-content-text"]/div[1]/table[6]'
demo=htmltab::htmltab(LINK,PATH)
demo=demo[,c(3,7:11)] # just what is needed
names(demo)=c("Country","elect","func","parti","cult","civil")
demo[,-1]=lapply(demo[,-1], as.numeric) 
demo=na.omit(demo)
h1=formula(elect~func + parti)
r1=lm(h1,data=demo)
summary(r1)

Call:
lm(formula = h1, data = demo)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7820 -1.2590 -0.0798  1.4786  4.4657 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.42114    0.43500  -5.566 1.04e-07 ***
func         0.73121    0.08212   8.904 9.73e-16 ***
parti        0.86349    0.10801   7.995 2.21e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.912 on 164 degrees of freedom
Multiple R-squared:  0.7531,    Adjusted R-squared:  0.7501 
F-statistic: 250.1 on 2 and 164 DF,  p-value: < 2.2e-16

2. El mismo investigador tiene una duda: cree que lograría un modelo mejor que el anterior si añadiese como predictor el nivel de libertades civiles. ¿Qué opinas?

h2=formula(elect~func + parti + civil)
r2=lm(h2,data=demo)

# model selection --------
anova(r1,r2)

Analysis of Variance Table

Model 1: elect ~ func + parti
Model 2: elect ~ func + parti + civil
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    164 599.32                                  
2    163 364.19  1    235.13 105.24 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

3. Descargue la data de este LINK (note que el documento incluye la metadata). No olvide trabajar con datos completos. ¿Se verifica la hipótesis que la Mortalidad materna es explicada por la fertilidad y la desnutrición?

rm(list = ls())
wb=read.csv("https://docs.google.com/spreadsheets/d/e/2PACX-1vSpe4ZRgc92WdHLjNdFJGDpEZOTZ2zvwVUUjPuRhZo-pXyoF44B91ut2ChagKB8Iw/pub?gid=1446884572&single=true&output=csv")
# names(wb) # fair
# str(wb)
wb=na.omit(wb)
h1=formula(MaterMort100m~TasaFertil1mMuje + UndernourishPerPop)
r1=lm(h1,data=wb)
summary(r1)

Call:
lm(formula = h1, data = wb)

Residuals:
    Min      1Q  Median      3Q     Max 
-254.49  -56.87  -20.80   16.02  941.71 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -12.1276    10.0810  -1.203  0.23055    
TasaFertil1mMuje     2.6156     0.2217  11.799  < 2e-16 ***
UndernourishPerPop   3.9128     1.0441   3.747  0.00024 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 134.7 on 180 degrees of freedom
Multiple R-squared:  0.9979,    Adjusted R-squared:  0.9979 
F-statistic: 4.335e+04 on 2 and 180 DF,  p-value: < 2.2e-16

# standardized results -------
lm.beta::lm.beta(r1)

Call:
lm(formula = h1, data = wb)

Standardized Coefficients::
       (Intercept)   TasaFertil1mMuje UndernourishPerPop 
                NA          0.7583632          0.2408552 

4. De lo anterior, Uds. se preguntan si pueden obtener un mejor modelo si añade el empleo, ¿Qué encuentras al correr el modelo actualizado?

h2=formula(MaterMort100m~TasaFertil1mMuje + UndernourishPerPop + EmployPerPop)
r2=lm(h2,data=wb)
anova(r1,r2)

Analysis of Variance Table

Model 1: MaterMort100m ~ TasaFertil1mMuje + UndernourishPerPop
Model 2: MaterMort100m ~ TasaFertil1mMuje + UndernourishPerPop + EmployPerPop
  Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
1    180 3267451                                  
2    179 2922157  1    345295 21.151 8.004e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

5. Diríjase a la pagina web de la CIA. Descargue Infant Mortality rate y Education expenditures (People and Society); Public debt (Economy); Carbon dioxide emissions (Energy). Ud quiere plantearse dos hipótesis. En ambas hipótesis, la deuda publica y las emisiones de carbono son variables independientes. En la primera hipótesis, la mortalidad infantil es la dependiente, y los gastos en educación pasa a ser independiente con las otras; en la segunda, los gastos en educación será la dependiente y la mortalidad infantil pasa a ser independiente con las ya mencionadas. Si Ud prueba ambas hipótesis, ¿Qué encuentras?

inf=htmltab::htmltab(infLink,xPathAll)

No encoding supplied: defaulting to UTF-8.
Warning: Columns [Date of Information] seem to have no data and are removed. Use rm_nodata_cols = F to suppress this behavior.

# hipotheses -------
h1=formula(inf~edu + deb + car)
h2=formula(edu~ inf + deb + car)
# modelling --------------------------------------------------------------
r1=lm(h1,data=alldata)
r2=lm(h2,data=alldata)
# results -------
summary(r1)
summary(r2)

6. Diversas escuelas de la ciudad están sufriendo quejas por ataques violentos a los alumnos. Los datos están disponibles en este ENLACE. Al probar la hipotesis que el tipo de escuela y la ubicación del colegio afectan los ataques ¿Qué encuentras?

rm(list = ls())
lk="https://docs.google.com/spreadsheets/d/e/2PACX-1vRZgLtdckN--CBTFFzWx_CZF8aqSosuwTZ6GOHAbBOZFEbIjnSvNg_0wrZ-VdW9T2u1qWdhogQIluYN/pub?gid=0&single=true&output=csv"
dataate=read.csv(lk)
#names(dataate)
modeltr <- glm(ataVio ~ tipo + region, family = poisson,
               offset = log(poblacion), data = dataate)
#summary(modeltr)
#overdispersion
AER::dispersiontest(modeltr,alternative='greater')$ p.value<0.05

[1] TRUE

#underdispersion
AER::dispersiontest(modeltr,alternative='less')$ p.value<0.05

[1] FALSE

library(MASS)
modeltr2 <- glm.nb(ataVio ~ tipo + region + offset(log(poblacion)), data = dataate)
summary(modeltr2)

Call:
glm.nb(formula = ataVio ~ tipo + region + offset(log(poblacion)), 
    data = dataate, init.theta = 1.242581365, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2499  -1.1136  -0.4544   0.4696   2.4795  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -8.20727    0.27203 -30.171  < 2e-16 ***
tipoPublica -0.16068    0.26537  -0.605  0.54486    
regionMW     0.12305    0.42403   0.290  0.77167    
regionNE     0.99469    0.34807   2.858  0.00427 ** 
regionS      0.54792    0.34520   1.587  0.11245    
regionW      0.06314    0.45582   0.139  0.88982    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(1.2426) family taken to be 1)

    Null deviance: 97.868  on 79  degrees of freedom
Residual deviance: 87.600  on 74  degrees of freedom
AIC: 420.59

Number of Fisher Scoring iterations: 1

              Theta:  1.243 
          Std. Err.:  0.276 

 2 x log-likelihood:  -406.588 

exp(coef(modeltr2))

 (Intercept)  tipoPublica     regionMW     regionNE      regionS      regionW 
0.0002726629 0.8515661558 1.1309419625 2.7038764497 1.7296441419 1.0651801125 

7. Se está estudiando los factores que influencian la admisión a los colegios elite del Perú. Para ello se tienen datos de los postulantes procedentes de diferentes escuelas primarias. Descargue la data de este LINK Si usa todas las variables disponibles como predictoras de la admisión, ¿Qué encuentras?

rm(list = ls())
lk="https://docs.google.com/spreadsheets/d/e/2PACX-1vSIeBRLEmHacvpWqAsi5W6B091eeZa0XN6wA3zww1qZlltK-EWCIkLN6J3_FkYvUDElViR67Rkx_80a/pub?gid=0&single=true&output=csv"
mydata <- read.csv(lk)
#str(mydata)
mylogit <- glm(admitido ~ letras + ciencias + prestigio, 
               data = mydata, family = "binomial")
#summary(mylogit)

sdVIs=apply(mydata[,c("letras","ciencias", "prestigio")],2,sd)
list(LogitSt=sdVIs*coef(mylogit)[c(2,3,4)])%>%
    data.frame()

summary(mylogit)

Call:
glm(formula = admitido ~ letras + ciencias + prestigio, family = "binomial", 
    data = mydata)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5802  -0.8848  -0.6382   1.1575   2.1732  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -6.249705   1.155732  -5.408 6.39e-08 ***
letras       0.002294   0.001092   2.101   0.0356 *  
ciencias     0.007770   0.003275   2.373   0.0177 *  
prestigio    0.560031   0.127137   4.405 1.06e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 499.98  on 399  degrees of freedom
Residual deviance: 459.44  on 396  degrees of freedom
AIC: 467.44

Number of Fisher Scoring iterations: 4

8. Se está estudiando qué factores influyen en que los convictos puestos en libertad vuelvan a la carcel. Los datos puede descargarlos de aqui. Con esos datos, dos investigadores desean probar algunas hipotesis. Un investigador “A” afirma que mientras mayor la edad, y si uno ya está casado, disminuirá el riesgo de volver a la carcel; mientras que otro investigador “B” asegura que, además de esas predictoras, la experiencia laboral previa también disminuye el riesgo de volver a la carcel. ¿Qué opinas en esta situación?

rm(list = ls())
lk="https://docs.google.com/spreadsheets/d/e/2PACX-1vQSlGaMI8Q8qlXI0Bp3m7BQcEh8ZLzaP7RymVtRYkg3ah1sZVlCi6-HmeKCic1RjfuH3gL_wrbMms88/pub?gid=1573532387&single=true&output=csv"
carcel=read.csv(lk)
#names(carcel)
library(survival)
carcel$survival=with(carcel,Surv(time = semanasLibre,event =  as.numeric(fueArrestado)))
carcel$casado=factor(carcel$casado)
COX_H1= formula(survival~edad + casado)
COX_H2= formula(survival~edad + casado + expLaboralPrevia)
#regression
rcox1 <- coxph(COX_H1,data=carcel)
rcox2 <- coxph(COX_H2,data=carcel)
summary(rcox1)

Call:
coxph(formula = COX_H1, data = carcel)

  n= 432, number of events= 114 

            coef exp(coef) se(coef)      z Pr(>|z|)   
edad    -0.06672   0.93545  0.02082 -3.205  0.00135 **
casado1 -0.49314   0.61070  0.37269 -1.323  0.18576   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

        exp(coef) exp(-coef) lower .95 upper .95
edad       0.9355      1.069    0.8981    0.9744
casado1    0.6107      1.637    0.2942    1.2678

Concordance= 0.616  (se = 0.027 )
Likelihood ratio test= 17.27  on 2 df,   p=2e-04
Wald test            = 13.99  on 2 df,   p=9e-04
Score (logrank) test = 14.51  on 2 df,   p=7e-04

summary(rcox2)

Call:
coxph(formula = COX_H2, data = carcel)

  n= 432, number of events= 114 

                     coef exp(coef) se(coef)      z Pr(>|z|)   
edad             -0.05614   0.94541  0.02158 -2.602  0.00927 **
casado1          -0.39569   0.67322  0.37911 -1.044  0.29662   
expLaboralPrevia -0.30907   0.73413  0.20427 -1.513  0.13026   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                 exp(coef) exp(-coef) lower .95 upper .95
edad                0.9454      1.058    0.9063    0.9862
casado1             0.6732      1.485    0.3202    1.4153
expLaboralPrevia    0.7341      1.362    0.4919    1.0956

Concordance= 0.615  (se = 0.026 )
Likelihood ratio test= 19.58  on 3 df,   p=2e-04
Wald test            = 16.83  on 3 df,   p=8e-04
Score (logrank) test = 17.64  on 3 df,   p=5e-04

anova(rcox1,rcox2)

Analysis of Deviance Table
 Cox model: response is  survival
 Model 1: ~ edad + casado
 Model 2: ~ edad + casado + expLaboralPrevia
   loglik  Chisq Df Pr(>|Chi|)
1 -666.75                     
2 -665.59 2.3132  1     0.1283

9. A su data de democracia añádale el índice de desarrollo humano (en ingles) que descargará de este LINK. En este caso quédese con los nombres que se escriban de la misma manera. Con la data unida, ¿Cuántos factores serían sugeridos de la data?

#-----
rm(list = ls())
library(htmltab)

# links
WhereDEMO=list(page="https://en.wikipedia.org/wiki/Democracy_Index",
               xpath='//*[@id="mw-content-text"]/div[1]/table[6]/tbody')

#carga
idh  = rio::import("HDR21-22_Statistical_Annex_HDI_Table.xlsx",skip=4)

New names:
• `` -> `...1`
• `` -> `...2`
• `` -> `...4`
• `` -> `...6`
• `` -> `...8`
• `` -> `...10`
• `` -> `...12`
• `` -> `...14`

demo  = htmltab(doc = WhereDEMO$page, 
                which  = WhereDEMO$xpath,
                encoding = "UTF-8")
#-----prepareIDH
#names(idh)
idh=idh[,c(2,3,5,7,9,11)]
newNames=c('Pais','puntuacion','EsperanzaVida','EscolaridadDuracion','EscolaridadPromedio','PBI')
names(idh)=newNames
#idh[idh$Pais=='Somalia' & !is.na(idh$Pais),]
idh=idh[c(1:202),]
idh=idh[!is.na(idh$Pais),]

idh[,-1]=lapply(idh[,-1], as.numeric)

Warning in lapply(idh[, -1], as.numeric) : NAs introduced by coercion
Warning in lapply(idh[, -1], as.numeric) : NAs introduced by coercion
Warning in lapply(idh[, -1], as.numeric) : NAs introduced by coercion
Warning in lapply(idh[, -1], as.numeric) : NAs introduced by coercion
Warning in lapply(idh[, -1], as.numeric) : NAs introduced by coercion

#idh[!complete.cases(idh[,-1]),]
idh=idh[complete.cases(idh[,-1]),]
row.names(idh)=NULL

#---- prepro DEMO
#names(demo)
demo=demo[,-c(1,2,6)]

names(demo)=c("Pais","RegimeType","Score","Electoral","Functioning","participation","culture",'Civilliberties')

demo[,3:8]=lapply(demo[,3:8],as.numeric)

#---- keys
idh$Pais= trimws(idh$Pais,whitespace = "[\\h\\v]")
demo$Pais= trimws(demo$Pais,whitespace = "[\\h\\v]") 

idhdemo=merge(idh,demo)

#----- efa

dontselect=c("Pais","RegimeType","puntuacion",'Score')
select=setdiff(names(idhdemo),dontselect) 
theData=idhdemo[,select]
library(polycor)
corMatrix=polycor::hetcor(theData)$correlations
library(psych)
#psych::KMO(corMatrix) 
#cortest.bartlett(corMatrix,n=nrow(theData))$p.value>0.05
library(matrixcalc)
#is.singular.matrix(corMatrix)
fa.parallel(theData, fa = 'fa',correct = T,plot = F)

Parallel analysis suggests that the number of factors =  2  and the number of components =  NA 

10. Si siguiese la recomendación para reducir la dimensionalidad, ¿Cuánta info logran recuperar los factores?

library(GPArotation)
resfa <- fa(theData,
            nfactors = 2,
            cor = 'mixed',
            rotate = "varimax",
            fm="minres")
print(resfa$loadings,cutoff = 0.5)

Loadings:
                    MR2   MR1  
EsperanzaVida       0.867      
EscolaridadDuracion 0.813      
EscolaridadPromedio 0.791      
PBI                 0.806      
Electoral                 0.901
Functioning         0.521 0.738
participation             0.767
culture                        
Civilliberties            0.913

                 MR2   MR1
SS loadings    3.450 3.443
Proportion Var 0.383 0.383
Cumulative Var 0.383 0.766

11. ¿Cómo interpreta el rol de complejidad y comunalidad de las variables utilizadas?

sort(resfa$communality)

            culture       participation EscolaridadPromedio 
          0.4264771           0.7059212           0.7288980 
                PBI         Functioning EscolaridadDuracion 
          0.7300781           0.8168660           0.8237711 
      EsperanzaVida           Electoral      Civilliberties 
          0.8384541           0.8777419           0.9443247 

sort(resfa$complexity)

          Electoral       EsperanzaVida                 PBI 
           1.160066            1.226160            1.242377 
     Civilliberties EscolaridadPromedio       participation 
           1.261717            1.318609            1.385843 
EscolaridadDuracion         Functioning             culture 
           1.463350            1.798043            1.983687 

12. Si plantease un análisis factorial basado en alguna teoría, ¿Qué usaría para darle sustento a los resultados de esa técnica?

modelCFA='demo=~ Electoral + Functioning + participation + culture + Civilliberties
idh=~EsperanzaVida + EscolaridadDuracion + EscolaridadPromedio + PBI'

# normalizar las variables:
theDataNorm=scale(theData)

library(lavaan)
cfa_fit <- cfa(modelCFA, data=theDataNorm, 
               std.lv=TRUE,  
               missing="fiml")
allParamCFA=parameterEstimates(cfa_fit,standardized = T)
allFitCFA=as.list(fitMeasures(cfa_fit))
allFitCFA[c("chisq", "df", "pvalue")] # pvalue>0.05

$chisq
[1] 120.3633

$df
[1] 26

$pvalue
[1] 4.274359e-14

allFitCFA$tli 

[1] 0.9017975

allFitCFA[c('rmsea.ci.lower','rmsea' ,'rmsea.ci.upper')] 

$rmsea.ci.lower
[1] 0.1302663

$rmsea
[1] 0.1582089

$rmsea.ci.upper
[1] 0.1873278

13. Con la data original del Democracy Index pide 4 clusters y compara el resultado con lo propuesto por The Economist. ¿Qué observas?

#### 1
dataClus=demo[,-c(1:3)]
row.names(dataClus)=demo$Pais
library(cluster)
g.dist = daisy(dataClus, metric="gower")
set.seed(123)
res.pam=pam(g.dist,4,cluster.only = F)

#nueva columna
dataClus$pam=res.pam$cluster

#aggregate(.~ pam, data=dataClus,mean)
demo$pam=res.pam$cluster
table(demo$pam,demo$RegimeType,dnn = c('Particion','Econo'))

         Econo
Particion Authoritarian Flawed democracy Full democracy Hybrid regime
        1             0               20             21             0
        2             0               25              0             6
        3             0                8              0            28
        4            59                0              0             0

14. Con la data original del IDH averigua cuantos clusters se recomiendan siguiendo la técnica aglomerativa.

#boxplot(idh[,-c(1:2)],horizontal = F,las=2,cex.axis = 0.5)


dataClus=BBmisc::normalize(idh[,-c(1:2)],method='standardize')
row.names(dataClus)=idh$Pais
library(cluster)
g.dist = daisy(dataClus, metric="gower")
set.seed(123)


library(factoextra)
fviz_nbclust(dataClus, hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F,hc_func = "agnes")

15. Con la data original del Democracy Index pide los clusters que se recomienden, ¿Cuántos países quedan mal clusterizados?

demo$pam=NULL
dataClus=demo[,-c(1:3)]
row.names(dataClus)=demo$Pais
library(cluster)
g.dist = daisy(dataClus, metric="gower")
fviz_nbclust(dataClus, pam,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F)

set.seed(123)
res.pam=pam(g.dist,7,cluster.only = F)
silPAM=data.frame(res.pam$silinfo$widths)
silPAM$country=row.names(silPAM)

poorPAM=silPAM[silPAM$sil_width<0,'country']#%>%sort()
poorPAM

 [1] "Trinidad and Tobago" "Slovakia"            "Kenya"              
 [4] "Tunisia"             "Ukraine"             "Albania"            
 [7] "Mexico"              "Guyana"              "Namibia"            
[10] "Ghana"               "Turkey"              "Hong Kong"          
[13] "Qatar"               "Angola"              "Eswatini"           
[16] "Laos"               

16. Haz un merge con IDH y Democracy Index, pero recuperando todos los paises aun cuando no se escriban igual. Normalice los valores del IDH con el metodo rango (del 1 al 10). ¿Hay diferencias en cuanto a las sugerencia de la cantidad de clusters para cada tecnica jerarquica si usa los componentes de ambos conceptos?

fviz_nbclust(dataClus, hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F,hc_func = "diana")

17. Para el caso anterior, si Ud pide 3 clusters, ¿Qué diferencias encuentra al clusterizar por ambas técnicas?

res.agnes<- hcut(g.dist, k = 3,hc_func='agnes',hc_method = "ward.D")
fviz_silhouette(res.agnes,print.summary = F)

res.diana <- hcut(g.dist, k = 3,hc_func='diana')
fviz_silhouette(res.diana,print.summary = F)

silDIANA=data.frame(res.diana$silinfo$widths)
silDIANA$country=row.names(silDIANA)
poorDIANA=silDIANA[silDIANA$sil_width<0,'country']#%>%sort()
poorDIANA

 [1] "Senegal"          "Liberia"          "Guatemala"       
 [4] "Turkey"           "Madagascar"       "Zambia"          
 [7] "Morocco"          "Malawi"           "Papua New Guinea"
[10] "Bolivia"          "Honduras"        

silAGNES=data.frame(res.agnes$silinfo$widths)
silAGNES$country=row.names(silAGNES)
poorAGNES=silAGNES[silAGNES$sil_width<0,'country']#%>%sort()
poorAGNES

 [1] "Guatemala"        "Madagascar"       "Turkey"          
 [4] "Zambia"           "Morocco"          "Malawi"          
 [7] "Papua New Guinea" "Honduras"         "Bolivia"         
[10] "Lesotho"          "Bhutan"           "Latvia"          
[13] "Lithuania"       

18. Rehaga las preguntas 16 y 17, pero normalizando toda la data con la técnica de estandarización.

idhdemo2=merge(idh,demo)
idhdemo2[,c(3:6,9:13)]=BBmisc::normalize(idhdemo2[,c(3:6,9:13)],method='standardize')
dataClus=idhdemo2[,c(3:6,9:13)]
row.names(dataClus)=idhdemo2$Pais

library(cluster)
g.dist = daisy(dataClus, metric="gower")
fviz_nbclust(dataClus, hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F,hc_func = "agnes")

fviz_nbclust(dataClus, hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F,hc_func = "diana")

res.agnes<- hcut(g.dist, k = 3,hc_func='agnes',hc_method = "ward.D")
fviz_silhouette(res.agnes,print.summary = F)

res.diana <- hcut(g.dist, k = 3,hc_func='diana')
fviz_silhouette(res.diana,print.summary = F)

silDIANA=data.frame(res.diana$silinfo$widths)
silDIANA$country=row.names(silDIANA)
poorDIANA=silDIANA[silDIANA$sil_width<0,'country']#%>%sort()
poorDIANA

 [1] "Senegal"          "Liberia"          "Guatemala"       
 [4] "Turkey"           "Madagascar"       "Zambia"          
 [7] "Morocco"          "Malawi"           "Papua New Guinea"
[10] "Bolivia"          "Honduras"        

silAGNES=data.frame(res.agnes$silinfo$widths)
silAGNES$country=row.names(silAGNES)
poorAGNES=silAGNES[silAGNES$sil_width<0,'country']#%>%sort()
poorAGNES

 [1] "Guatemala"        "Madagascar"       "Turkey"          
 [4] "Zambia"           "Morocco"          "Malawi"          
 [7] "Papua New Guinea" "Honduras"         "Bolivia"         
[10] "Lesotho"          "Bhutan"           "Latvia"          
[13] "Lithuania"