#+++++++++++++++++++++++++++++++++++++++++
# Laboratorio 1- ANALISI DI DATI SU RETE -  Marzo 2025
#+++++++++++++++++++++++++++++++++++++++++

## Dati: Esportazioni tra paesi 


# libraries ---------------------------------------------------------------
install.packages("amen")
library(amen)


# dataset -----------------------------------------------------------------
data(IR90s)
?IR90s
str(IR90s)
attributes(IR90s$dyadvars)$dimnames[[3]]  # variabili diadiche
IR90s$dyadvars[1:10,1:10,2]               # non simmetrica
IR90s$dyadvars[1:10,1:10,3]               # simmetrica

attributes(IR90s$nodevars)$dimnames[[2]]  # variabili di nodo
IR90s$nodevars[1:10,1] #considero la prima
IR90s$nodevars #posso guardarle tutte, noto che polity assume valori positivi e negativi

esport <- IR90s$dyadvars[,,2] #creo una variabile 'esport'

library(lattice)
levelplot(esport) #rappresentazione in forma di matrice

esport[1:10,1:10]
levelplot(esport[1:10,1:10])


# dataset -----------------------------------------------------------------
#data(IR90s)
#Y sono le connessioni dove c'e' almeno una esportazione
# se è vera la proposizione >0 allora attribuiamo il valore 1 

Y<-(esport>0)*1 #moltiplico per 1 per trasfomare in numerica una matrice che altrimenti sarebbe 'logica' (vero/falso)
levelplot(Y)


#seleziono i primi dieci paesi in ordine alfabetico
Y10<- Y[1:10, 1:10]
levelplot(Y10) 
#possiamo renderci conto facilmente di quali paesi sono maggiori esportatori e importatori (osserviamo ad esempio il Belgio) 

#----
#selezione solo dei 30 paesi con gdp maggiore
gdp<-IR90s$nodevars[,2]
sort(gdp, decreasing=T)
topgdp<-which(gdp>=sort(gdp,decreasing=TRUE)[30])
Y30<-Y[topgdp,topgdp] #matrice dei top gdp

levelplot(Y30) #è interessante notare come qui tutte le relazioni sono positive

##############voglio sapere se c'è una connessione in entrata o in uscita
###Creazione di matrice di adiacenza indiretta
#creo matrice di adiacenza indiretta
#considero i due stati connessi
Yu<- ((Y+t(Y))>0) #matrice solo di vero/falso
Yu<- ((Y+t(Y))>0)[1:10,1:10]
Yu<- ((Y+t(Y))>0)+0 
#anche qui faccio una operazione per trasformare la matrice: aggiungo lo zero per trasfomare in numerica una matrice che altrimenti sarebbe 'logica' (vero/falso)

levelplot(Yu)

#
isSymmetric(Yu)
isSymmetric(Y)
isSymmetric(Y30)



#carico pacchetto per grafici e statistiche di rete
install.packages("igraph")
library(igraph)
##creo dei grafi che poi uso per calcolare statistiche di rete e grafici 

net = graph_from_adjacency_matrix(Y, mode = 'directed', weighted = NULL, diag=F)
#str(net)

net10 = graph_from_adjacency_matrix(Y10, mode = 'directed', weighted = NULL, diag=F)

net30 = graph_from_adjacency_matrix(Y30, mode = 'directed', weighted = NULL, diag=F)


netu <-  graph_from_adjacency_matrix(Yu, mode = 'undirected', weighted = NULL, diag=F)
#str(netu)




### Analisi dei dati in quanto rete
####tutta la rete
E(net)
V(net)

edge_density(net)


#####rete dei 30 paesi con GDP più elevato
E(net30)
V(net30)

edge_density(net30)


#####rete indiretta
E(netu)
V(netu)

edge_density(netu)



## distribuzione del grado
degree(net)
hist(degree(net),breaks = 10,col='lavender')
which.max(degree(net))

degree(net30)
hist(degree(net30),col='lavender')


degree(netu)
hist(degree(netu),breaks = 10,col='lavender')


#?shortest_paths
shortest_paths(net, from = 'ITA')
shortest_paths(net, from = 'USA')

shortest_paths(net, from = 'ITA', to='AFG')
shortest_paths(net, from = 'USA', to='CUB')
shortest_paths(net, from = 'UGA', to='CUB')

#######rete dei 30 paesi 
shortest_paths(net30, from = 'COL', to='EGY')
shortest_paths(net30, from = 'IRN', to='EGY')
shortest_paths(net30, from = 'EGY', to='IRN')

shortest_paths(netu, from = 'EGY', to='IRN')




#altre misure viste a lezione
hist(betweenness(net), breaks = 20, col = 'lavender')

hist(sqrt(betweenness(net)), breaks = 20, col = 'lavender')
diameter(net)

set.seed(111) 
#se voglio che il grafico resti sempre uguale, altrimenti se lo rifaccio ogni volta cambia configurazione. Posso anche provare varie configurazioni, e trovare quella che più mi soddisfa. 

plot(net)

#per rete dei 30
#
hist(betweenness(net30), breaks = 20, col = 'lavender')

hist(sqrt(betweenness(net30)), breaks = 20, col = 'lavender')
diameter(net30)

plot(net30)

plot(net10)
#per rete indiretta
hist(betweenness(netu), breaks = 20, col = 'lavender')

hist(sqrt(betweenness(netu)), breaks = 20, col = 'lavender')
diameter(netu)

plot(netu)

#######################################################################
### Analisi con  modelli per dati di rete- Marzo 2025
###Dyadic Data analysis with amen (https://arxiv.org/abs/1506.08237)
# libraries ---------------------------------------------------------------
#install.packages("amen")
library(amen)


# dataset -----------------------------------------------------------------
data(IR90s)
str(IR90s)


#### ---- dati di esportazione per i 30 paesi con maggiore GDP
gdp<-IR90s$nodevars[,2]
topgdp<-which(gdp>=sort(gdp,decreasing=TRUE)[30] )
 
Y<-log( IR90s$dyadvars[topgdp,topgdp,2] + 1)
#hist(Y)
Y[1:5,1:5]

# ANOVA -------------------------------------------------------------------
#### ---- ANOVA 
Rowcountry<-matrix(rownames(Y),nrow(Y),ncol(Y))
Rowcountry[1:5,1:5]

Colcountry<-t(Rowcountry)
Colcountry[1:5,1:5]

mod = lm( c(Y) ~ c(Rowcountry) + c(Colcountry), na.action = na.exclude) 
anova(mod)


#### ---- confronto tra paesi
rmean<-rowMeans(Y,na.rm=TRUE) 
cmean<-colMeans(Y,na.rm=TRUE)
muhat<-mean(Y,na.rm=TRUE)
ahat<-rmean-muhat
bhat<-cmean-muhat

# effetti di riga
ahat[1:5]
head( sort(ahat,decreasing=TRUE) )
# effetti di colonna
bhat[1:5]
head( sort(bhat,decreasing=TRUE) )

#### ---- matrice di varianze e covarianze
cov(cbind(ahat,bhat) ) 
cor(cbind(ahat,bhat) )
#cor( ahat, bhat)

# visualizzo correlazioni
plot(ahat, bhat, type="n")
text(ahat,bhat,names(ahat))
abline(h=0,lty=2,col="grey")
abline(v=0,lty=2,col="grey")
abline(a=0,b=1,lty=2,col="grey")




# Social Relation Model ---------------------------------------------------
fit_SRM<-ame(Y, family="nrm")
muhat              # media generale
mean(fit_SRM$BETA) # SRM-media generale
summary(fit_SRM)

fit_SRM$APM      #stima degli effetti di riga
fit_SRM$BPM      #stima degli effetti di colonna


gofstats(Y) 
     
# statistiche osservate (deviazione standard empirica di medie di riga, deviazione standard empirica di medie di colonna, correlazione diadica empirica, dipendenza triadica empirica)
#Gli istogrammi rappresentano i valori delle gofastats simulate e vanno confrontati con le statistiche osservate gofstats(Y) (righe verticali rosse)

plot(fit_SRM)



# Social Relation Regression Model ----------------------------------------
#### ---- variabili di nodo
dimnames(IR90s$nodevars)[[2]]
Xn<-IR90s$nodevars[topgdp,]
Xn[,1:2]<-log(Xn[,1:2])

#### ---- variabili diadiche
dimnames(IR90s$dyadvars)[[3]]
Xd<-IR90s$dyadvars[topgdp,topgdp,c(1,3,4,5)]   # tutte le variabili diadiche eccetto export
Xd[,,2]<-log(Xd[,,2])

#### ---- Modello
fit_SRRM<-ame(Y,Xd=Xd,Xr=Xn,Xc=Xn, family="nrm")
summary(fit_SRRM)

#### ---- Regressione lineare
fit_rm<-ame(Y,Xd=Xd,Xr=Xn,Xc=Xn,family="nrm", rvar=FALSE,cvar=FALSE,dcor=FALSE)
summary(fit_rm)



#### ---- Modello AME con 2 fattori latenti
fit_ame2<-ame(Y,Xd,Xn,Xn,R=2, family="nrm")
summary(fit_ame2)
circplot(Y,U=fit_ame2$U,V=fit_ame2$V)

#il circle plot può aiutare a identificare gruppi di nodi che sono simili in termini di esportazioni e importazioni.