load("Prima_elaborazione_3.RData")
library(igraph)
library(HistDAWass)
library(brainGraph)
library(tidyverse)
library(plotly)
library(effects)
library(heatmaply)
source("Preprocessing.R")
source("useful_functions.R")
knitr::opts_chunk$set(cache=TRUE)

PCA of 4500 nets

Plane 1-2

plo12<-ggplotly(plot_indiv(res_ALL,x=1,y=2)+geom_point(aes(color=as.factor(DESC$types))))
plo12

Plane 3-4

plo34<-ggplotly(plot_indiv(res_ALL,x=3,y=4)+geom_point(aes(color=as.factor(DESC$types))))
plo34              

a 3d plot

DF3d<-as.data.frame(res_ALL$ind$coord)
DF3d<-DF3d %>% mutate(ID=row.names(.),types=as.factor(DESC$types))
plo3d<- plot_ly(DF3d, x = ~Dim.1, y = ~Dim.2, z = ~Dim.3,color = ~types,
                text=~ID,
                type = 'scatter3d', mode = 'markers', 
                marker=list(
                  size=4,
                  line=list(
                    color='MediumPurple',
                    width=0.5
                  )
                ))
plo3d

Cluster analysis

tree=hclust(dist(res_ALL$ind$coord),method = "ward.D")
clu<-cutree(tree,45)

Compute the Gini index for each class vs clusters

DF<-data.frame(types=DESC$types,size=DESC$size,dens=DESC$dens)

DF<-DF %>% mutate(exper=paste(types,size,dens,sep="_"))
print(paste0("Adjusted Rand index -> ",mclust::adjustedRandIndex(clu,DF$exper)))
## [1] "Adjusted Rand index -> 0.677877853119289"
confu<-as.data.frame(table(DF$exper,clu))
confu %>% pivot_wider(names_from = clu,values_from = Freq)
confu<-confu %>% pivot_wider(names_from = clu,values_from = Freq)
pur_cla<-numeric();for(i in 2:46){pur_cla<-c(pur_cla,sum((confu[,i]/sum(confu[,i]))^2))}
pur_cla_row<-numeric();for(i in 1:45){
  pur_cla_row<-c(pur_cla_row,sum((confu[i,2:46]/sum(confu[i,2:46]))^2))
}

Effects

tmp<-data.frame(CLU=confu$Var1,pp=pur_cla_row)
tmp<-tmp %>% right_join(unique(DF), by=c("CLU"="exper"))
tmp$types<-factor(tmp$types, levels=c("RG","2_I","3_I","PA","SW"))
tmp$size<-factor(tmp$size, levels=c(960,480,240))
tmp$dens<-factor(tmp$dens, levels=c(0.1,0.05,0.01))

model<-lm(pp~types+size+dens,tmp)
print(summary(model))
## 
## Call:
## lm(formula = pp ~ types + size + dens, data = tmp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3147 -0.1136  0.0129  0.1025  0.3227 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.911542   0.077494  11.763 6.85e-14 ***
## types2_I    -0.251178   0.081686  -3.075 0.004004 ** 
## types3_I    -0.312378   0.081686  -3.824 0.000502 ***
## typesPA     -0.096867   0.081686  -1.186 0.243455    
## typesSW      0.166222   0.081686   2.035 0.049277 *  
## size480     -0.050200   0.063274  -0.793 0.432757    
## size240      0.008093   0.063274   0.128 0.898932    
## dens0.05    -0.010920   0.063274  -0.173 0.863946    
## dens0.01    -0.180267   0.063274  -2.849 0.007208 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1733 on 36 degrees of freedom
## Multiple R-squared:  0.6079, Adjusted R-squared:  0.5208 
## F-statistic: 6.977 on 8 and 36 DF,  p-value: 1.589e-05

Plots of effects

plot(allEffects(model))

Heatmap of clusters

confum<-as.matrix(confu[,2:46])

row.names(confum)<-confu$Var1
plo_hm<-heatmaply(confum,Colv=NA)
plo_hm