Introducción

Una red dinámica hace referencia a una red que cambia a lo largo del tiempo. Del instante \(t\) al \(t+1\) pueden aparecer o desaparecer nodos o aristas.

El adjetivo “dinámico” se ha aplicado en el contexto de redes para describir tres aspectos diferentes de un sistema en evolución:

  • Dinámicas de aristas: Formación y disolución de enlaces
  • Dinámicas de nodos: Entrada y salida de actores
  • Dinámicas en los atributos: Atributos de los nodos y las aristas pueden cambiar en el tiempo.

Si ambos tipos de dinámica se presentan juntas se dice que la red co-evoluciona.

Conceptualmente, definimos una red dinámica como un grafo indexado en el tiempo \[G(t)=\left(V(t), E(t)\right)\]

con \(t\) variando de forma discreta o continua.

  • \(V(t)\) como el conjunto de vértices presente en la red en el momento \(t\).
  • \(E(t)\) como el conjunto de aristas \(e_{ij}(t)\), para \(i,j\in V(t)\) que indica presencia o ausencia de aristas entre los vértices.

En varias ocasiones es conveniente dejar el conjunto de vértices constante a lo largo del tiempo, en cuyo caso escribimos \(V(t)=V\) para todo \(t\).

Tipos de datos

La información de una red temporal se puede recopilar de diferentes maneras, lo que conduce a varios tipos diferentes de datos:

  • Datos tipo panel representando redes en tiempo discreto (redes longitudinales)
  • Datos transversales de la red con información sobre la duración del enlace
  • Redes temporales en tiempo continuo

Tipos de modelos

Hay tres marcos generales de modelado para el análisis de datos de redes temporales en el campo de las redes sociales.

  • Modelos estocásticos orientados al actor (SAOM): permiten modelar la co-evolución de atributos nodales y dinámicas de aristas. Además, se formulan como modelos de tiempo continuo, aunque típicamente se estiman par datos de tiempo discreto.

  • Modelos de eventos relacionales: son modelos de tiempo continuo para la dinámica de los nodos y los enlaces, por lo que pueden asumir la independencia de las díadas y explotar la metodología de la regresión logística.

  • Modelos de grafos aleatorios exponenciales temporales (TERGMs): representan dinámicas de los enlaces únicamente, y se formulan en tiempo discreto.

Representación y visualización

Con los diferentes tipos de redes dinámicas mencionadas vienen así mismo varias formas en las que estas redes pueden visualizarse.

Panel

Bajo esta representación es importante mantener cierta “estabilidad” a lo largo de las visualizaciones manteniendo las posiciones de los vértices invariantes.

Consideremos los 3 paneles de los datos Sampson.

library(networkDynamic)
library(igraph)
library(ergm)  
data(samplk)

samp.list <- list(samplk1,samplk2,samplk3)
samp.dyn <- networkDynamic(network.list = samp.list, start = 1)
## Onsets and termini not specified, assuming each network in network.list should have a discrete spell of length 1
## Argument base.net not specified, using first element of network.list instead
## Created net.obs.period to describe network
##  Network observation period info:
##   Number of observation spells: 1 
##   Maximal time range observed: 1 until 4 
##   Temporal mode: discrete 
##   Time unit: step 
##   Suggested time increment: 1
library(sna)
coord <- gplot.layout.fruchtermanreingold(
  as.matrix.network.adjacency(samp.dyn), layout.par = list())

par(mfrow = c(2,2), oma=c(1,1,1,1), mar=c(4,1,3,1))

plot(network.extract(samp.dyn, at = 1), coord = coord,main = "Tiempo 1", 
     displaylabels = T, label.cex = 0.6, vertex.cex = 2, pad = 0)
plot(network.extract(samp.dyn, at = 2), coord = coord,main = "Tiempo 2", 
     displaylabels = T, label.cex = 0.6, vertex.cex = 2, pad = 0)
plot(network.extract(samp.dyn, at = 3), coord = coord,main = "Tiempo 3", 
     displaylabels = T, label.cex = 0.6, vertex.cex = 2, pad = 0)
plot(samp.dyn, coord = coord,main = "Colapsada", 
     displaylabels = T, label.cex = 0.6, vertex.cex = 2, pad = 0)

Esta representación, aunque sencilla, no es la más útil si queremos información sobre los cambios temporales en las amistades entre monjes del monasterio, así como cuando los datos son lo suficientemente ricos temporalmente que graficar cada ‘slice’ se vuelve insostenible.

Película

Para obtener un rápido resumen visual de cómo cambia la estructura con el tiempo, podemos renderizar la red dinámica como una animación

library(ndtv)
# Datos incluidos en el paquete
data(short.stergm.sim)

render.d3movie(short.stergm.sim,
               plot.par=list(displaylabels=T),
               output.mode = 'htmlWidget')
## slice parameters:
##   start:0
##   end:25
##   interval:1
##   aggregate.dur:1
##   rule:latest

Este método, no obstante, no se recomienda del todo para redes muy granes, debido al costo computacional.

Línea de tiempo aproximada

Otra herramienta interesante de visualización es el cronograma de proximidad. En este gráfico, los nodos se posicionan en el eje-y de acuerdo a su proximidad en términos de distancias geodésicas. Es útil ya que los cambios de posición de los nodos nos habla de cambios en la estructura de la red.

proximity.timeline(short.stergm.sim,default.dist = 6,
          mode = 'sammon',labels.at = 17,vertex.cex = 4)
## collapsing slice networks ...
## computing vertex positions using 1D sammon layout ...
##   computing positions for slice 1
##   computing positions for slice 2
##   computing positions for slice 3
##   computing positions for slice 4
##   computing positions for slice 5
##   computing positions for slice 6
##   computing positions for slice 7
##   computing positions for slice 8
##   computing positions for slice 9
##   computing positions for slice 10
##   computing positions for slice 11
##   computing positions for slice 12
##   computing positions for slice 13
##   computing positions for slice 14
##   computing positions for slice 15
##   computing positions for slice 16
##   computing positions for slice 17
##   computing positions for slice 18
##   computing positions for slice 19
##   computing positions for slice 20
##   computing positions for slice 21
##   computing positions for slice 22
##   computing positions for slice 23
##   computing positions for slice 24
##   computing positions for slice 25
## rendering splines for each vertex ...

Note que es fácil identificar que la familia ‘Lamberteschi’ se desprende de su grupo de a partir de \(t=15\), volviéndose un nodo aislado por aproximadamente 5 periodos de tiempo.

Nota: El eje \(y\) de este gráfico toma valores negativos ya que este representa la posición en un espacio unideimensional que aproxima las distancias geodésicas de la red. Nos cercanos en \(y\) tienden a estar próximos en la red, y el signo (+ o -) no tiene interpretación, ya que lo que importa son las diferencias entre posiciones. Por ejemplo, un conjunto de nodos A, B y C toma los valores en el eje \(y\) de -4, -3 y 5. Esto se interpreta como que A y B están cerca geodésicamente y C está alejado de ambos.

Análisis descriptivo

Veamos un ejemplo de análisis descriptivo de una red dinámica nos basaremos en Kolaczyk (2020).

El objeto hc contiene información de 75 personas acerca de registros de contactos cara a cara entre pacientes y diversos tipos de trabajadores de la salud en la unidad geriátrica de un hospital en Lyon, Francia, en 2010, desde la 1 p. m. del lunes 6 de diciembre hasta las 2 p. m. del viernes 10 de diciembre.

Un objetivo principal de este estudio fue obtener información sobre el patrón de contactos en un entorno hospitalario de este tipo y la forma en que podría transmitirse la infección.

# Datos
library(sand)
data(hc)
head(hc)
##   Time ID1 ID2  S1  S2
## 1  140  15  31 MED ADM
## 2  160  15  22 MED MED
## 3  500  15  16 MED MED
## 4  520  15  16 MED MED
## 5  560  16  22 MED MED
## 6  580  16  22 MED MED
# Numero de personas por status
ID.stack <- c(hc$ID1,hc$ID2)
Status.stack <- c(as.character(hc$S1), as.character(hc$S2))
my.t <- table(ID.stack,Status.stack)
v.status <- character(nrow(my.t))

for(i in (1:length(v.status))){
  v.status[i] <- names(which(my.t[i,]!=0))}

table(v.status)
## v.status
## ADM MED NUR PAT 
##   8  11  27  29
# Numero de relaciones entre status
status.t <- table(hc$S1,hc$S2)
status.t <- status.t + t(status.t)
diag(status.t) <- round(diag(status.t)/2)

status.t
##      
##         ADM   MED   NUR   PAT
##   ADM   279   459  2596   441
##   MED   459  5660  1769  1471
##   NUR  2596  1769 12695  6845
##   PAT   441  1471  6845   209

Representación:

tmp.es <- paste(hc$S1, "-", hc$S2, sep = "")

e.status <- character(dim(hc)[[1]])
e.status[tmp.es=="ADM-ADM"] <- "ADM-ADM"
e.status[tmp.es=="MED-MED"] <- "MED-MED"
e.status[tmp.es=="NUR-NUR"] <- "NUR-NUR"
e.status[tmp.es=="PAT-PAT"] <- "PAT-PAT"
e.status[(tmp.es=="ADM-MED") | (tmp.es=="MED-ADM")] <- "ADM-MED"
e.status[(tmp.es=="ADM-NUR") | (tmp.es=="NUR-ADM")] <- "ADM-NUR"
e.status[(tmp.es=="ADM-PAT") | (tmp.es=="PAT-ADM")] <- "ADM-PAT"
e.status[(tmp.es=="MED-NUR") | (tmp.es=="NUR-MED")] <- "MED-NUR"
e.status[(tmp.es=="MED-PAT") | (tmp.es=="PAT-MED")] <- "MED-PAT"
e.status[(tmp.es=="NUR-PAT") | (tmp.es=="PAT-NUR")] <- "NUR-PAT"

my.hc <- data.frame(Time = hc$Time/(60*60), ID1 = hc$ID1, ID2 = hc$ID2,
                    Status = e.status)

# Histogramas
library(lattice)
histogram(~Time|Status, data = my.hc, xlab="Horas", layout = c(5,2),
          ylab = "Porcentaje del total")

Histogramas de contactos en el conjunto de datos hospitalarios hc, agrupados según el estado del par que interactúa.

Distribución del grado

Veamos la caracterización de una red dinámica utilizando los datos de contacto del hospital, donde, el objetivo es comprender cómo podría transmitirse la infección entre individuos en este entorno.

# Distribucion del grado por cada periodo de tiempo

all.deg <- sapply(g.sl12, igraph::degree)
sl.lab <- sapply(1:8, function(i) paste(12*(i-1), "-", 12*i, "hrs", sep=""))
deg.df <- data.frame(Degree = as.vector(all.deg), 
                     Slice  = rep(sl.lab,each=75),
                     Estado = rep(V(g.week)$Status, times=8))

library(ggplot2)
p = qplot(factor(Degree), data = deg.df, geom = "bar", fill = Estado)
p + facet_grid(Slice~.) + xlab("Grado") + ylab("Frec.")

Distribuciones de grados para los datos de contacto del hospital, una por cada una de las ocho redes estáticas correspondientes a periodos consecutivos de 12 horas

  • Se aprecia una moda importante en cero en cada gráfico de barras. Esto es, durante cualquier período de 12 horas, una fracción considerable de las 75 personas estudiadas no participó en ningún contacto registrado.

  • La gran mayoría de los pacientes tienen contactos con un número comparativamente menor de personas.

  • El personal del hospital tiende a tener contactos con un mayor número de personas.

# Por grado

# Top 5 de nodos de mayor grado por periodo
top.deg <- lapply(1:8,function(i) {all.deg[,i][rank(all.deg[,i])>=70]})

# Numero de periodos en los que un nodo estuvo en el top 5 de mayor grado
tab <- table(unlist(lapply(1:8,function(i)
  as.numeric(names(top.deg[[i]])))))

# Nodos y numero de periodos de 12h en el top 5 de mayor grado
sort(tab, decreasing = TRUE)
## 
##  7 15 17 19 23  1  5 11 13 21 24 26 27 29 37 64  8 10 22 25 31 34 36 63 
##  4  4  3  3  3  2  2  2  2  2  2  2  2  2  2  2  1  1  1  1  1  1  1  1

Los nodos 7 y 15 estuvieron en el top 5 de mayor grado en cuatro de los ocho periodos consecutivos de 12 horas. Como se verá a continuación, estos nodos hacen referencia a una enfermera y un médico, respectivamente.

# Identificar nodos
V(g.week)$Status[c(7,15)]
## [1] "NUR" "MED"
# Tiempo promedio por contacto en cada periodo de 12h
all.str <- sapply(g.sl12, strength)

# Tiempo de contacto/ grado
all.r <- all.str/all.deg
Tiempo promedio por contacto en cada periodo
Nodo Periodo 1 Periodo 2 Periodo 3 Periodo 4 Periodo 5 Periodo 6 Periodo 7 Periodo 8
7 2.00 25.79 11.1 32.73 14.20 33.19 8.33 37.34
15 29.71 26.33 17.0 12.48 19.27 15.30 19.40 12.93 
# Comparación
summary(c(all.r))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    1.00    5.00   10.43   12.06   15.73   47.35     261

Esta enfermera y este médico tienden a tener una tasa de contactos igual o superior a la media, es decir, entran en contacto con muchas más personas de lo habitual en esta sala del hospital.

En consecuencia, estos dos individuos son candidatos naturales para ser considerados como potenciales propagadores.

# Distancia promedio a los demas nodos por cada periodo de tiempo
sp.len <- lapply(1:8, function(i) {
  spl <- distances(g.sl12[[i]],v = c(7,15), to = V(g.sl12[[i]]), weights = NA)
  spl[spl==Inf] <- NA
  spl
  })

# Tabla
ave.spl <- sapply(1:8,function(i) apply(sp.len[[i]], 1, mean, na.rm=T))
Distancia promedio a los demás nodos por cada periodo de tiempo
Nodo Periodo 1 Periodo 2 Periodo 3 Periodo 4 Periodo 5 Periodo 6 Periodo 7 Periodo 8
7 3.05 1.27 1.79 1.12 1.8 1.33 2 1.24
15 1.72 1.51 2.48 1.35 1.36 1.26 1.61 1.36 
# Distancia promedio de los caminos mas cortos por periodo de tiempo
round(sapply(g.sl12,mean_distance),digits=2)
## [1] 3.80 3.26 4.47 4.33 3.69 3.56 4.06 3.50

La enfermera y el médico parecen estar más cerca de los demás en la red la mayor parte del tiempo, lo que podría tener un efecto perjudicial en la contención de la infección, si alguno de ellos se contagiara.

Modelos TERGM

Los modelos TERGMs son extensiones de los ERGMs para el modelamiento de redes dinámicas en tiempo discreto.

Los TERGMs postulan dos o más modelos para las dinámicas de las aristas. Cada uno de estos es, por si solo, un ERGM que opera sobre un conjunto específico de nodos y enlaces.

Modelo conjunto vs Modelo separable

Al abordar los Modelos Temporales de Grafos Aleatorios Exponenciales (TERGMs), el punto de partida fundamental es decidir cómo conceptualizamos el cambio de la red de un momento a otro. Naturalmente, existen dos marcos de modelado principales que definen esta transición: los Modelos Conjuntos y los Modelos Separables (STERGMs)

Modelo conjunto

Los modelos conjuntos representan la evolución de la red como un proceso único y unificado. En este enfoque, se utiliza una sola fórmula para explicar el paso de la red del tiempo \(t\) al tiempo \(t+1\). La suposición subyacente es que los mismos procesos sociales y factores estructurales (como la reciprocidad o la transitividad) impulsan tanto la creación de nuevos vínculos como la terminación de los existentes.

STERGM

“Separable” aquí significa que se supone que la formación de aristas es independiente de la disolución de estas dentro de cada paso del tiempo (\(t\) a \(t+1\)), y dependiente de Markov entre los pasos (la formación de una arista en \(t\) puede influenciar a la disolución de una arista adyacente en \(t+1\)). Esto permite que los factores que influyen en la formación sean diferentes de aquellos que influyen en la disolución.

En el contexto de las amistades, esto significa que, estadísticamente, el modelo asume que el proceso de “quiénes se hacen amigos hoy” no tiene nada que ver con el proceso de los que están “peleando y rompiendo hoy”. Así mismo, si bien la formación y la disolución son independientes dentro de un paso temporal, hay dependencia entre los pasos temporales.

Formalmente considere la evolución de una red en el tiempo \(t\) al tiempo \(t+1\).

Se define la red de formación \(\mathbf{Y}^+\) como aquella que contiene todos los enlaces existentes en \(\mathbf{Y}^t\) junto con los que se forman entre \(t\) y \(t+1\). Así: \[\mathbf{Y}^t\subseteq \mathbf{Y}^+,\] donde \(\mathbf{Y}^t\) es la red inicial.

Dada \(\mathbf{Y}^t = \boldsymbol{y}^t\), la red de formación \(\mathbf{Y}^+\) es generada de un ERGM:

\[\textsf{Pr}(\mathbf{Y}^+=\boldsymbol{y}^+|\mathbf{Y}^{t};\boldsymbol{\theta}^+)=\frac{\exp\{\boldsymbol{\theta}^{+\top}\mathbf{g}^+(\boldsymbol{y}^+)\}}{\kappa}, \quad \boldsymbol{y^+}\in \mathcal{Y}^+(\boldsymbol{y}^t)\] donde:

  • \(\boldsymbol{\theta}^+\) es el vector de parámetros de formación.
  • \(\mathbf{g}^+(\boldsymbol{y}^+)\) es el vector de estadísticas asociadas a la formación de la red.
  • \(\mathcal{Y}^+(\boldsymbol{y}^t)\) es el espacio de las redes obtenibles agregando enlaces a \(\boldsymbol{y}^t\).
  • \(\kappa\) es la constante normalizadora.

La red de disolución \(\mathbf{Y}^-\) (bajo un abuso de notación) se puede representar como \(\mathbf{Y}^-=\mathbf{Y}^t-D\), donde \(D\) es el conjunto de enlaces que se disuelven de la red entre \(t\) y \(t+1\): \[\mathbf{Y}^-\subseteq \mathbf{Y}^t.\]

Esta es simultáneamente generada por un ERGM, como sigue: \[\textsf{Pr}(\mathbf{Y}^-=\boldsymbol{y}^-|\mathbf{Y}^{t};\boldsymbol{\theta}^-)=\frac{\exp\{\boldsymbol{\theta}^{-\top}\mathbf{g}^-(\boldsymbol{y}^-)\}}{\kappa}, \quad \boldsymbol{y^-}\in \mathcal{Y}^-(\boldsymbol{y}^t),\] donde:

  • \(\boldsymbol{\theta}^-\) es el vector de parámetros de disolución
  • \(\mathbf{g}^-(\boldsymbol{y}^-)\) es el vector de estadísticas asociadas a la disolución de la red.
  • \(\mathcal{Y}^-(\boldsymbol{y}^t)\) es el espacio de las redes obtenibles eliminando enlaces de \(\boldsymbol{y}^t\).

Nota: \(\mathbf{Y}^-\) contiene los enlaces que sobreviven tras la disolución, no los enlaces que desaparecen.

La red observada en el siguiente periodo se obtiene agregando los nuevos enlaces formados y eliminando los enlaces disueltos: \[\mathbf{Y}^{t+1}=\mathbf{Y}^t \cup (\mathbf{Y}^+-\mathbf{Y}^t)-(\mathbf{Y}^t-\mathbf{Y}^-). \]

Dado que \(\mathbf{Y}^t\subseteq\mathbf{Y}^+\), la expresión anterior puede simplificarse a:

\[\mathbf{Y}^{t+1}=\mathbf{Y}^+-(\mathbf{Y}^t-\mathbf{Y^-}), \] o equivalentemente: \[\mathbf{Y}^{t+1}=\mathbf{Y}^-\cup(\mathbf{Y}^+-\mathbf{Y}^t), \]

Si \(\mathbf{Y}^+\) es condicionalmente independiente de \(\mathbf{Y}^-\) dada \(\mathbf{Y}^{t}\) entonces la probabilidad de transición se expresa como: \[\textsf{Pr}(\mathbf{Y}^{t+1}=\boldsymbol{y}^{t+1}|\mathbf{Y}^t=\boldsymbol{y}^t)= \textsf{Pr}(\mathbf{Y}^+=\boldsymbol{y}^+|\mathbf{Y}^t=\boldsymbol{y}^t; \boldsymbol{\theta}^+)\times \textsf{Pr}(\mathbf{Y}^-=\boldsymbol{y}^-|\mathbf{Y}^t=\boldsymbol{y}^t; \boldsymbol{\theta}^-)\] lo que permite estimar separadamente los parámetros de formación y de disolución.

Operadores temporales

Una transición entre dos redes: \[\mathbf{Y}^{t}\longrightarrow\mathbf{Y}^{t+1},\] puede analizarse desde varias perspectivas:

  • ¿Qué enlaces nuevos aparecieron?
  • ¿Qué enlaces desaparecieron?
  • ¿Qué proporción de la red final está ocupada?
  • ¿Cuánto cambió la red en total?

Los operadores Form(), Persist(), Diss(), Change() y Cross() responden a estas preguntas, representando un tipo específico de dinámica de los enlaces.

Operador A qué responde
Form() De todas las díadas que podían formar un enlace, ¿cuáles efectivamente lo formaron?
Persist() De los vínculos que existen en \(t\), ¿cuáles lograron sobrevivir hasta \(t+1\)?
Diss() De los enlaces existentes, ¿cuáles desaparecieron en \(t+1\)?
Change() ¿Qué tan probable es que un par de nodos cambie su estado entre dos momentos?
Cross() ¿Qué enlaces existen en \(t+1\)?

Form()

Modela exclusivamente la formación de enlaces entre tiempos. Es por esto que solo considera los pares de nodos que están vacíos en \(t\), ignorando los enlaces ya existentes. Esta se puede interpretar como \[\textsf{Pr}(Y_{ij}^{t+1}=1|Y_{ij}^{t}=0)=\frac{\text{enlaces formados}}{\text{díadas vacías en }t}\]

Persist()

Modela la permanencia de los enlaces existentes de \(t\) a \(t+1\). Esta se interpreta como \[\textsf{Pr}(Y_{ij}^{t+1}=1|Y_{ij}^{t}=1)=\frac{\text{enlaces persistentes}}{\text{enlaces existentes en }t}\]

Diss()

Modela exclusivamente la disolución de enlaces entre tiempos. Razón por la cual solo considera los pares de nodos no vacíos en \(t\). Se interpreta como: \[\textsf{Pr}(Y_{ij}^{t+1}=0|Y_{ij}^{t}=1)=\frac{\text{enlaces disueltos}}{\text{enlaces existentes en }t}\]

Este operador es complementario con persist(), ya que \(\text{Pr(persistencia)}+ \text{Pr(disolución)}=1\). Así mismo, un coeficiente positivo indica que cierta característica incrementa la probabilidad de ruptura del enlace.

Change()

Mide la inestabilidad general de la red. No importa si el cambio fue una formación o una disolución; sólo importa que el estado del vínculo haya cambiado. Se interpreta como \[\textsf{Pr}(Y_{ij}^{t+1}\neq Y_{ij}^{t})=\frac{\text{díadas que cambiaron}}{\text{todas las díadas}}\]

Cross()

Este último operador es el más diferente de todos, ya que no modela transiciones sino a la red observada directamente en el tiempo \(t+1\). Este es simplemente la densidad de la red en el tiempo \(t+1\), en su caso más simple: \[\textsf{Pr}(Y_{ij}^{t+1}=1)=\frac{\text{enlaces en }t+1}{\text{todas las díadas posibles}}\]

Ignora completamente el estado previo de la red y se concentra únicamente en explicar la estructura observada en el instante actual, de ahí viene su nombre cross-sectional (transversal).

Si bien los operadores se pueden mezclar de diferentes maneras, hay dos combinaciones naturales que se utilizan con mayor frecuencia:

  • Form() + Diss() (un modelo ‘separable’)

  • Cross() + Change() (un modelo ‘conjunto’)

Ejemplos

Hay dos opciones para el argumento estimate en un modeo TERGM:

  • CMLE (estimación condicional de máxima verosimilitud): para ajustar un modelo a redes dinámicas tipo panel

  • EGMME (método generalizado de estimación de momentos en equilibrio): apara ajustar un modelo a una sola red transversal con información de la duración de los enlaces.

CMLE

Para datos tipo panel

La forma canónica de datos para modelar procesos dinámicos de red consiste en observaciones de una red completa en dos o más puntos en el tiempo en el mismo conjunto de nodos.

Como primer ejemplo, se muestra el modelo separable Form() + Persist() (este último se elije sobre el operador de disolución ya que el signo de persistencia tiene el mismo impacto que el operador de formación)

library(tergm)


set.seed(123)
samp.fit.fp <- tergm(samp.list ~
    Form(~edges+mutual+cyclicalties+transitiveties) +
    Persist(~edges+mutual+cyclicalties+transitiveties),
  estimate = "CMLE",
  times = c(1:3)
    )

summary(samp.fit.fp)
## Call:
## tergm(formula = samp.list ~ Form(~edges + mutual + cyclicalties + 
##     transitiveties) + Persist(~edges + mutual + cyclicalties + 
##     transitiveties), estimate = "CMLE", times = c(1:3))
## 
## Monte Carlo Conditional Maximum Likelihood Results:
## 
##                           Estimate Std. Error MCMC % z value Pr(>|z|)    
## Form(1)~edges              -3.4766     0.3450      0 -10.076   <1e-04 ***
## Form(1)~mutual              2.0370     0.4160      0   4.896   <1e-04 ***
## Form(1)~cyclicalties       -0.1366     0.1949      0  -0.701   0.4833    
## Form(1)~transitiveties      0.3915     0.2391      0   1.637   0.1016    
## Persist(1)~edges            0.2035     0.3055      0   0.666   0.5054    
## Persist(1)~mutual           0.8017     0.5214      0   1.537   0.1242    
## Persist(1)~cyclicalties    -0.2036     0.2542      0  -0.801   0.4232    
## Persist(1)~transitiveties   0.5309     0.2647      0   2.005   0.0449 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 848.4  on 612  degrees of freedom
##  Residual Deviance: 375.4  on 604  degrees of freedom
##  
## AIC: 391.4  BIC: 426.7  (Smaller is better. MC Std. Err. = 1.537)

No hay muchos términos significativos aquí. Esto se debe en parte a que esta es una pequeña red dispersa, con solo dos cambios observados, por lo que la cantidad de información disponible para probar estos efectos es mínima.

Vemos que una relación es más propensa a formarse si en \(t+1\) cierra una díada mutua. El log-odds condicional incrementa en un 2.04, lo que quiere decir que la probabilidad relativa de que se forme un vínculo es 7.67 mayor si ese vínculo va a completar una relación recíproca, en comparación con uno que no lo hace.

Los coeficientes para tríadas transitivas son positivos y levemente significativos para formación y persistencia. Los coeficientes para tríadas cíclicas son negativos, pero no significativos.

set.seed(123)
samp.fit.cc <- tergm(samp.list ~
    Cross(~edges+mutual+cyclicalties+transitiveties) +
    Change(~edges+mutual+cyclicalties+transitiveties),
  estimate = "CMLE",
  times = c(1:3)
    )

summary(samp.fit.cc)
## Call:
## tergm(formula = samp.list ~ Cross(~edges + mutual + cyclicalties + 
##     transitiveties) + Change(~edges + mutual + cyclicalties + 
##     transitiveties), estimate = "CMLE", times = c(1:3))
## 
## Monte Carlo Conditional Maximum Likelihood Results:
## 
##                          Estimate Std. Error MCMC % z value Pr(>|z|)    
## Cross(1)~edges           -1.89693    0.24203      0  -7.837   <1e-04 ***
## Cross(1)~mutual           1.65992    0.39289      0   4.225   <1e-04 ***
## Cross(1)~cyclicalties    -0.03647    0.17644      0  -0.207   0.8362    
## Cross(1)~transitiveties   0.46345    0.20483      0   2.263   0.0237 *  
## Change(1)~edges          -1.65568    0.23339      0  -7.094   <1e-04 ***
## Change(1)~mutual         -0.28897    0.50759      0  -0.569   0.5691    
## Change(1)~cyclicalties   -0.09086    0.18241      0  -0.498   0.6184    
## Change(1)~transitiveties  0.46223    0.24422      0   1.893   0.0584 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 848.4  on 612  degrees of freedom
##  Residual Deviance: 379.5  on 604  degrees of freedom
##  
## AIC: 395.5  BIC: 430.8  (Smaller is better. MC Std. Err. = 1.394)

Nuevamente, una arista es más propensa a estar presente en la red si es parte de una díada mutua; el log-odds condicional aumenta un 1.66 lo que se traduce en que es 5.36 veces más probable encontrar un vínculo si este es mutuo. Recordemos que Cross describe la estrutura de la red, por lo que es corecto pensar que hay mucha reciprocidad en la red, sin importar el intervalo de tiempo en la que se mire. No sorprende esta interpretación ya que en el primer modelo las díadas mutuas tienden a formarse y a persistir en el tiempo.

El efecto de la reciprocidad en Change es negativo aunque no significativo. Si fuese significativo se interpretaría como que los vínculos sociales que forman parte de un par mutuo son más estables en el tiempo, que los que no lo son.

EGMME

Para redes transversales

Ahora, imaginamos otro escenario en el que se observó una sola red transversal con información retrospectiva sobre la duración de las aristas. ¿Cómo estimar los modelos de formación y persistencia cuando solo se cuenta con una red? … Imponiendo supuestos de equilibrio:

  • Balance de flujos: la tasa a la que se forman nuevos vínculos es igual a la tasa a la que desaparecen los enlaces existentes (disolución), por lo que ambos procesos se compensan entre si.

  • Estabilidad en las estadísticas: las propiedades globales de la red (densidad total, número de triángulos, reciprocidad, etc.) permanecen constantes a lo largo del tiempo

Si se asume que un proceso está en equilibrio, se tiene que \[\text{prevalencia}=\text{incidencia}\times\text{duración},\]

donde la prevalencia es la cantidad total de aristas existentes en la red transversal, la incidencia se refiere a la tasa de aparición de nuevos vínculos, y la duración es la información que se tiene sobre las aristas de ante mano.

La incidencia es lo que alimenta la red. Si conoces la prevalencia (la foto) y la duración de los vínculos (la persistencia), el modelo EGMME puede calcular matemáticamente cuál tuvo que ser la incidencia (tasa de formación) necesaria para mantener la red en ese nivel de equilibrio

Ejemplo: Red de amistad

  1. Considere que se estudia una red de amistad y, por encuestas, se sabe que una amistad dura, en promedio, 10 años (donde cada paso de tiempo \(t\) es de 1 año).

    En un modelo simple para la disolución, solo usando edges, existe una relación directa entre la probabilidad de persistencia \(p\) y la duración esperada \(d\): \(d=\frac{1}{1-p}\)

    Si \(d=10\), entonces: \[10=\frac{1}{1-p}\rightarrow 1-p= 0.1 \rightarrow p=0.9\] Esto significa que, para que una amistad dure 10 años, cada año debe tener un \(90\%\) de probabilidad de sobrevivir.

    Este se convierte a escala log-odds para que sea el coeficiente \(\theta_{\text{pers}}\) que conocemos: \[\theta_{\text{persistencia}}=\ln \left(\frac{p}{1-p}\right)=\ln\left(\frac{0.9}{0.1}\right)=\ln(9)\approx 2.197\]

    Nota: El coeficiente de disolución no es nada más que \(-2.197\)

  2. Ya sabiendo que tan rápido “mueren” las amistades (disolución), el modelo debe encontrar ahora que tan rápido se forman para que la red se mantenga igual a la ‘foto’ observada.

    Supongase que la red transversal tiene una densidad del \(5\%\) (prevalencia). El modelo EGMME se pregunta: “Si estoy perdiendo el 10% de mis vínculos en cada paso de tiempo, ¿cuántos vínculos nuevos debo crear para que la densidad se quede estancada en el 5%?”.

    El modelo, empieza con coeficientes de formación al azar, con estos simula redes y calcula sus densidades. Sigue este proceso, hasta acercarse lo más posible al \(5\%\).

Modelos SAOM

SAOM: Modelo Estocástico Orientado al Actor. (Luke, 2015)

  • Se basa en un proceso de Markov, es decir, la probabilidad de que una red cambie de estado depende únicamente de su estado actual.

  • Los enlaces pueden variar.

  • Queremos distinguir efectos de la influencia social y la selección social.

Datos

Los datos de Coevolve se presentan en forma de una lista de cuatro redes igraph.

Se trata de redes de amistad entre 37 estudiantes y estado de tabaquismo en cuatro momentos del tiempo.

Análisis descriptivo

library(igraph)
library(UserNetR)
library(RSiena)

# Datos de red de amistad (4 periodos de tiempo)
data(Coevolve)

fr_w1 <- igraph::upgrade_graph(Coevolve$fr_w1)
fr_w2 <- igraph::upgrade_graph(Coevolve$fr_w2)
fr_w3 <- igraph::upgrade_graph(Coevolve$fr_w3)
fr_w4 <- igraph::upgrade_graph(Coevolve$fr_w4)
# Ordenes
c(vcount(fr_w1), vcount(fr_w2), vcount(fr_w3), vcount(fr_w4))
## [1] 37 37 37 37
# Tamaños 
c(ecount(fr_w1), ecount(fr_w2), ecount(fr_w3), ecount(fr_w4))
## [1] 179 179 179 179

Visualización:

# Grafico

colors <- c("darkgreen", "lightskyblue2") # azul: fumador, verde: no fumador
shapes <- c("circle", "square")           # circulo: mujer, cuadrado: hombre

# Layout Kamada-Kawai
coord  <- layout_with_kk(fr_w1)

# Visualizacion (4 periodos de tiempo)

par(mfrow = c(2,2), mar = c(1,1,2,9), xpd = TRUE)

plot(fr_w1, vertex.color = colors[V(fr_w1)$smoke+1],
     vertex.shape = shapes[V(fr_w1)$gender],
     vertex.size = 10, main = "Oleada 1", vertex.label = NA,
     edge.arrow.size = 0.2, layout = coord)

plot(fr_w2, vertex.color = colors[V(fr_w2)$smoke+1],
     vertex.shape = shapes[V(fr_w2)$gender],
     vertex.size = 10, main = "Oleada 2", vertex.label = NA,
     edge.arrow.size = 0.2, layout = coord)

plot(fr_w3, vertex.color = colors[V(fr_w3)$smoke+1],
     vertex.shape = shapes[V(fr_w3)$gender],
     vertex.size = 10, main = "Oleada 3", vertex.label = NA,
     edge.arrow.size = 0.2, layout = coord)

plot(fr_w4, vertex.color = colors[V(fr_w4)$smoke+1],
     vertex.shape = shapes[V(fr_w4)$gender],
     vertex.size = 10, main = "Oleada 4", vertex.label = NA,
     edge.arrow.size = 0.2, layout = coord)

# Leyenda lateral
legend("right", inset = c(-0.45, 0),
       legend = c("No fumador", "Fumador", "Mujer", "Hombre"),
       col = c("darkgreen", "lightskyblue2", "black", "black"),
       pch = c(16, 16, 1, 0), pt.cex = 1.5, bty = "n")

Cambios en el tiempo de la red Coevolve (Índice de Jaccard: 0.8)

Podemos observar que la estructura topológica de la red es similar cuando pasamos del tiempo \(t\) a \(t+1\), pues el índice de Jaccard, el cual mide la similitud entre redes, es de 0.8, como no es exactamente 1, las redes no son idénticas. Este es un escenario adecuado para ajuste del modelo SAOM, pues el cambio es gradual, si estamos en presencia de cambios abruptos (Jaccard menor a 0.3), diremos que el cambio no es gradual, y esto dificultará el análisis.

Características de la red Coevolve a lo largo de cuatro oleadas
Oleada Orden Densidad Grado.ent..medio Fumadores Modularidad
1 37 0.134 4.838 8 0.001
2 37 0.134 4.838 10 0.044
3 37 0.134 4.838 12 0.077
4 37 0.134 4.838 14 0.129

Podemos que el estado de fumador aumenta con el tiempo de forma lineal. Es más, sigue exactamente este patrón:

\[ \text{Fumadores}(t) = 2t + 6, \quad t = 1,2,3,4 \]

Preparar datos

# Convertir matrices de adyacencia a tipo "matrix"
matw1 <- as.matrix(as_adjacency_matrix(fr_w1))
matw2 <- as.matrix(as_adjacency_matrix(fr_w2))
matw3 <- as.matrix(as_adjacency_matrix(fr_w3))
matw4 <- as.matrix(as_adjacency_matrix(fr_w4))
# Son las mismas matrices?
identical(matw1 , matw2)
## [1] FALSE

Veamos la especificación de las variables previo al ajuste del modelo SAOM:

  • Variable de comportamiento (dependiente): atributo de un nodo que cambia con el tiempo (estado de fumador)

  • coCovar: atributo de un nodo que no cambia con el tiempo (género).

  • varCovar: atributo que sí cambia con el tiempo.

  • Dependiente: lazos de amistad.

Nota: las matrices de adyacencia deben ser tipo matrix cuando se construya el arreglo.

# Apilar matrices

# n = 37 individuos
# 4 momentos en el tiempo

# Arreglo con matrices de adyacencia
fr4wav <- sienaDependent(array(c(matw1, matw2, matw3, matw4),
                               dim = c(37,37,4)), sparse = FALSE)

# Genero (constante)
gender_vect <- V(fr_w1)$gender
gender      <- coCovar(gender_vect, centered = FALSE)

# Estado de fumador (comportamiento)
smoke <- array(c(V(fr_w1)$smoke,V(fr_w2)$smoke,
                 V(fr_w3)$smoke,V(fr_w4)$smoke), dim = c(37,4))

smokebeh <- sienaDependent(smoke, type = "behavior")

# Lazos de amistad (dependiente)
friend <- sienaDataCreate(fr4wav, smokebeh, gender)
# Verificar cambio gradual
# Indice Jaccard > 0.3 (En este caso J = 0.817)

# print01Report(friend, modelname = 'Coevolve Example') (genera .txt)

Especificación del modelo

Dado que estamos explorando un modelo coevolutivo, examinaremos los efectos sobre la probabilidad de formación de vínculos (fr4wav) y los efectos sobre los cambios en el comportamiento (smokebeh).

# Especificacion de efectos
frndeff <- getEffects(friend)
# Ver effectsDocumentation(frndeff) (824 efectos)
# HIPOTESIS: (6 en total)

# 1. Homofilia de genero
frndeff <- includeEffects(frndeff, sameX,
                          interaction1 = "gender", name = "fr4wav")
##   effectNumber effectName  shortName include fix   test  initialValue parm
## 1 336          same gender sameX     TRUE    FALSE FALSE          0   0
# 2. Según el hábito de fumar del individuo, este tiene mayor o menor probabilidad 
# de formar una amistad (egoX). (seleccion social) 
frndeff <- includeEffects(frndeff, egoX,
                          interaction1 = "smokebeh", name = "fr4wav")
##   effectNumber effectName   shortName include fix   test  initialValue parm
## 1 464          smokebeh ego egoX      TRUE    FALSE FALSE          0   0
# 3. La probabilidad de formar una amistad esta relacionada con el hábito de fumar de los compañeros (altX).
frndeff <- includeEffects(frndeff, altX,
                          interaction1 = "smokebeh", name = "fr4wav")
##   effectNumber effectName     shortName include fix   test  initialValue parm
## 1 449          smokebeh alter altX      TRUE    FALSE FALSE          0   0
# 4. Homofilia de tabaquismo
frndeff <- includeEffects(frndeff, sameX,
                          interaction1 = "smokebeh", name = "fr4wav")
##   effectNumber effectName    shortName include fix   test  initialValue parm
## 1 520          same smokebeh sameX     TRUE    FALSE FALSE          0   0
# 5. La probabilidad de cambiar de comportamiento está relacionada con la similitud 
# promedio del hábito de fumar entre todos los individuos vinculados (avSim) (influencia social)
frndeff <- includeEffects(frndeff, avSim,
                          interaction1 = "fr4wav", name = "smokebeh")
##   effectNumber effectName                  shortName include fix   test 
## 1 724          smokebeh average similarity avSim     TRUE    FALSE FALSE
##   initialValue parm
## 1          0   0
# 6. Tendencia a que los lazos de amistad sean recíprocos
frndeff <- includeEffects(frndeff, recip,
                          name = "fr4wav")
##   effectNumber effectName  shortName include fix   test  initialValue parm
## 1 18           reciprocity recip     TRUE    FALSE FALSE          0   0

Estimación del modelo

# Estimacion del modelo
myalgorithm <- sienaAlgorithmCreate(projname = 'coevolve')
## If you use this algorithm object, siena07 will create/use an output file coevolve.txt .
# Semilla
set.seed(999)

# nbrNodes = 3: usa 3 de 4 nucleos del computador
# returnDeps = TRUE: para ver graficos de bondad de ajuste y simulaciones

RSmod1 <- siena07(myalgorithm, data = friend,
                  effects = frndeff, batch = TRUE,
                  verbose = FALSE, useCluster = TRUE,
                  initC = TRUE, nbrNodes = 3, returnDeps = TRUE)
## 
## Start phase 0 
## theta:  2.004  2.004  2.004 -0.808  0.000  0.000  0.000  0.000  0.000  0.208  0.208  0.208  0.562  0.000 
## 
## Start phase 1 
## Phase 1 Iteration 1 Progress: 0%
## Phase 1 Iteration 4 Progress: 0%
## Phase 1 Iteration 10 Progress: 0%
## Phase 1 Iteration 25 Progress: 0%
## Phase 1 Iteration 40 Progress: 0%
## theta:  1.7875  1.8964  1.8750 -0.8904  0.1693  0.0482  0.0449 -0.0250  0.0070  0.2393  0.2190  0.2155  0.6579  1.0000 
## 
## Start phase 2.1
## Phase 2 Subphase 1 Iteration 1 Progress: 25%
## Phase 2 Subphase 1 Iteration 2 Progress: 25%
## theta  1.4302  1.7547  1.5825 -1.0773  0.4704  0.1674  0.1111 -0.0806  0.0420  0.2724  0.2410  0.2104  0.8269  1.9288 
## ac  2.283  0.930  1.908  1.184  1.960  1.634 -3.128  1.039  1.260 -0.600 -5.859 -0.478 -1.138  1.329 
## Phase 2 Subphase 1 Iteration 3 Progress: 25%
## Phase 2 Subphase 1 Iteration 4 Progress: 25%
## theta  0.878  1.341  1.114 -1.523  0.955  0.456  0.277 -0.195  0.188  0.350  0.276  0.266  1.682  4.482 
## ac  1.8568  1.0195  1.8172  1.1326  1.9105  1.6650 -2.9685  1.0488  1.2742  0.0311 -1.7520 -0.3367 -1.0310  1.0752 
## Phase 2 Subphase 1 Iteration 5 Progress: 25%
## Phase 2 Subphase 1 Iteration 6 Progress: 25%
## theta  0.782  1.131  1.071 -1.894  1.184  0.675  0.406 -0.235  0.359  0.366  0.316  0.340  2.484  6.889 
## ac  1.1546  1.0328  1.5656  1.1522  1.5150  1.6760 -2.2616  1.0487  0.9603 -0.0842 -1.3836 -0.3329 -0.7847  1.1461 
## Phase 2 Subphase 1 Iteration 7 Progress: 25%
## Phase 2 Subphase 1 Iteration 8 Progress: 25%
## theta  0.944  1.189  1.145 -2.225  1.244  0.882  0.476 -0.271  0.539  0.410  0.334  0.365  3.121  8.672 
## ac  1.232  0.961  1.559  1.136  1.269  1.567 -0.951  1.005  1.009 -0.330 -1.304 -0.238 -0.225  0.754 
## Phase 2 Subphase 1 Iteration 9 Progress: 25%
## Phase 2 Subphase 1 Iteration 10 Progress: 26%
## theta  1.153  1.023  1.060 -2.463  1.070  1.010  0.537 -0.294  0.793  0.314  0.292  0.355  3.347  8.892 
## ac  1.151  0.955  1.537  1.125  1.236  1.505 -1.410  1.026  1.032 -0.503 -0.844 -0.196 -0.559  0.758 
## theta  1.1404  1.1221  1.1336 -3.1354  0.8443  1.4590  0.7464 -0.0793  1.2447  0.2997  0.3748  0.3609  6.9419 19.5822 
## ac  0.0915  0.1746  0.3693  0.7998  0.8163  0.4477  0.0370  0.0740  0.0519 -0.2711 -0.0967 -0.0396  0.0179  0.1147 
## theta:  1.1404  1.1221  1.1336 -3.1354  0.8443  1.4590  0.7464 -0.0793  1.2447  0.2997  0.3748  0.3609  6.9419 19.5822 
## 
## Start phase 2.2
## Phase 2 Subphase 2 Iteration 1 Progress: 33%
## Phase 2 Subphase 2 Iteration 2 Progress: 33%
## Phase 2 Subphase 2 Iteration 3 Progress: 33%
## Phase 2 Subphase 2 Iteration 4 Progress: 33%
## Phase 2 Subphase 2 Iteration 5 Progress: 33%
## Phase 2 Subphase 2 Iteration 6 Progress: 33%
## Phase 2 Subphase 2 Iteration 7 Progress: 33%
## Phase 2 Subphase 2 Iteration 8 Progress: 33%
## Phase 2 Subphase 2 Iteration 9 Progress: 33%
## Phase 2 Subphase 2 Iteration 10 Progress: 33%
## theta  1.1346  1.1320  1.1211 -3.0716  0.8380  1.4184  0.7165 -0.0711  1.2134  0.3101  0.3699  0.3480  7.1181 19.8062 
## ac -0.10100 -0.02267 -0.07416  0.00865  0.08278  0.04741 -0.02965  0.15427  0.02463 -0.17435  0.01673  0.01701 -0.01058  0.12173 
## theta:  1.1346  1.1320  1.1211 -3.0716  0.8380  1.4184  0.7165 -0.0711  1.2134  0.3101  0.3699  0.3480  7.1181 19.8062 
## 
## Start phase 2.3
## Phase 2 Subphase 3 Iteration 1 Progress: 41%
## Phase 2 Subphase 3 Iteration 2 Progress: 41%
## Phase 2 Subphase 3 Iteration 3 Progress: 41%
## Phase 2 Subphase 3 Iteration 4 Progress: 41%
## Phase 2 Subphase 3 Iteration 5 Progress: 41%
## Phase 2 Subphase 3 Iteration 6 Progress: 41%
## Phase 2 Subphase 3 Iteration 7 Progress: 41%
## Phase 2 Subphase 3 Iteration 8 Progress: 41%
## Phase 2 Subphase 3 Iteration 9 Progress: 41%
## Phase 2 Subphase 3 Iteration 10 Progress: 41%
## theta  1.1443  1.1373  1.1346 -3.0478  0.8421  1.4179  0.6979 -0.0856  1.1774  0.3027  0.3730  0.3497  7.4799 20.2595 
## ac  0.02697 -0.04932  0.00656  0.11967 -0.07553  0.12176 -0.00985  0.04249  0.00895  0.05614  0.03492 -0.02556 -0.01613  0.05810 
## theta:  1.1443  1.1373  1.1346 -3.0478  0.8421  1.4179  0.6979 -0.0856  1.1774  0.3027  0.3730  0.3497  7.4799 20.2595 
## 
## Start phase 2.4
## Phase 2 Subphase 4 Iteration 1 Progress: 51%
## Phase 2 Subphase 4 Iteration 2 Progress: 51%
## Phase 2 Subphase 4 Iteration 3 Progress: 51%
## Phase 2 Subphase 4 Iteration 4 Progress: 51%
## Phase 2 Subphase 4 Iteration 5 Progress: 51%
## Phase 2 Subphase 4 Iteration 6 Progress: 51%
## Phase 2 Subphase 4 Iteration 7 Progress: 51%
## Phase 2 Subphase 4 Iteration 8 Progress: 51%
## Phase 2 Subphase 4 Iteration 9 Progress: 51%
## Phase 2 Subphase 4 Iteration 10 Progress: 51%
## theta  1.1411  1.1344  1.1262 -3.0678  0.8368  1.4221  0.7149 -0.0754  1.2050  0.3108  0.3765  0.3561  6.9482 19.9796 
## ac -0.04466  0.05859  0.12891 -0.09815 -0.01110 -0.10218  0.02771  0.04015  0.04192  0.01481  0.04692 -0.00575 -0.04456 -0.04922 
## theta:  1.1411  1.1344  1.1262 -3.0678  0.8368  1.4221  0.7149 -0.0754  1.2050  0.3108  0.3765  0.3561  6.9482 19.9796 
## 
## Start phase 3 
## Phase 3 Iteration 1 Progress 67%

Cómo interpretar:

  • En la columna Estimate veremos cómo contribuye el efecto de interés a la probabilidad de formación de enlaces. Si el signo es positivo, la probabilidad aumenta, si es negativo, disminuye.

  • Los t-ratio dan cuenta de la falta de convergencia, es decir, esperamos que sean pequeños (menores a 0.15).

  • Para determinar si un efecto es significativo o no, tendremos que calcular \(\frac{\text{Estimate}}{\text{Std. Error}}\) el cual se distribuye \(t\), pero asintóticamente se comporta como una normal estándar, luego se verá si es mayor que 1.96 o no. Si lo es, diremos que el efecto es significativo al 95% de confianza.

# Resumen del modelo
summary(RSmod1)
## Estimates, standard errors and convergence t-ratios
## 
##                                            Estimate   Standard   Convergence 
##                                                         Error      t-ratio   
## Network Dynamics 
##    1. rate constant fr4wav rate (period 1)  1.1411  (   0.2009 )   -0.0080   
##    2. rate constant fr4wav rate (period 2)  1.1344  (   0.2589 )    0.0753   
##    3. rate constant fr4wav rate (period 3)  1.1262  (   0.2054 )    0.0197   
##    4. eval outdegree (density)             -3.0678  (   0.4092 )   -0.0073   
##    5. eval reciprocity                      0.8368  (   0.2365 )   -0.0118   
##    6. eval same gender                      1.4221  (   0.3029 )   -0.0192   
##    7. eval smokebeh alter                   0.7149  (   0.3748 )    0.0252   
##    8. eval smokebeh ego                    -0.0754  (   0.3003 )   -0.0017   
##    9. eval same smokebeh                    1.2050  (   0.4437 )   -0.0057   
## 
## Behavior Dynamics
##   10. rate rate smokebeh (period 1)         0.3108  (   0.1782 )   -0.0158   
##   11. rate rate smokebeh (period 2)         0.3765  (   0.2742 )    0.0306   
##   12. rate rate smokebeh (period 3)         0.3561  (   0.1936 )    0.0305   
##   13. eval smokebeh linear shape            6.9482  (  49.4802 )    0.0049   
##   14. eval smokebeh average similarity     19.9796  ( 149.7879 )   -0.0139   
## 
## Overall maximum convergence ratio:    0.1029 
## 
## 
## Total of 2281 iteration steps.
## 
## Covariance matrix of estimates (correlations below diagonal)
## 
##        0.040        0.005       -0.001        0.001        0.000        0.003       -0.002       -0.001       -0.004        0.003       -0.002       -0.003       -0.732       -2.287
##        0.101        0.067       -0.001        0.010       -0.003       -0.005        0.001       -0.011       -0.003        0.008       -0.031       -0.013       -7.523      -23.090
##       -0.031       -0.017        0.042        0.005       -0.001       -0.002       -0.003        0.001       -0.004        0.001        0.000       -0.004        0.172        0.351
##        0.006        0.095        0.065        0.167       -0.002       -0.079       -0.068       -0.020       -0.133        0.004       -0.014       -0.004       -1.116       -3.511
##        0.002       -0.042       -0.023       -0.018        0.056       -0.011       -0.004       -0.009       -0.012       -0.001        0.000       -0.001        0.039        0.197
##        0.043       -0.068       -0.037       -0.634       -0.150        0.092        0.000        0.005        0.006       -0.004        0.015        0.007        2.983        9.187
##       -0.022        0.006       -0.036       -0.446       -0.049       -0.003        0.141       -0.002        0.101        0.000       -0.012       -0.007       -2.771       -8.626
##       -0.022       -0.139        0.014       -0.164       -0.124        0.056       -0.018        0.090        0.032       -0.006        0.016        0.004        3.749       11.342
##       -0.042       -0.022       -0.039       -0.732       -0.113        0.042        0.608        0.244        0.197        0.002       -0.006       -0.002       -3.019       -9.279
##        0.073        0.175        0.016        0.057       -0.013       -0.078        0.004       -0.118        0.031        0.032       -0.010       -0.004       -2.922       -8.708
##       -0.043       -0.435       -0.005       -0.125        0.004        0.185       -0.113        0.193       -0.047       -0.212        0.075        0.015        8.570       26.457
##       -0.078       -0.258       -0.105       -0.051       -0.030        0.111       -0.093        0.077       -0.028       -0.115        0.292        0.037        3.819       11.832
##       -0.074       -0.587        0.017       -0.055        0.003        0.199       -0.149        0.252       -0.138       -0.331        0.632        0.399     2448.290     7399.333
##       -0.076       -0.595        0.011       -0.057        0.006        0.202       -0.154        0.252       -0.140       -0.326        0.644        0.408        0.998    22436.421
## 
## Derivative matrix of expected statistics X by parameters:
## 
##       29.221        0.000        0.000        1.510       -4.767       -1.360       -0.205       -1.430        4.353       -0.338        0.000        0.000       -0.447        0.502
##        0.000       28.603        0.000       -0.523       -3.546       -2.925        1.075       -1.716        2.786        0.000        0.042        0.000        0.217       -0.233
##        0.000        0.000       27.656       -0.110       -1.706       -0.620        1.527       -0.238        3.090        0.000        0.000       -0.475       -1.024        0.617
##        2.226        0.475        0.638      100.665       64.649       85.932       -2.646       -4.755       66.215        0.019        0.620        0.280       -0.072        0.913
##       -0.939       -0.458        0.434       32.722       59.599       30.843       -1.132       -0.236       21.909        0.264        0.228        0.267        0.975        0.327
##        0.529       -0.875        0.488       85.957       61.254       85.161       -2.372       -3.521       55.966        0.254        0.592        0.214       -0.179        0.980
##        0.077        0.661        0.501        1.085       -0.482        0.705       15.543        2.637       -7.869        0.308        0.456        0.135        1.039       -0.353
##        0.114       -0.200        0.674       -2.670        1.619       -1.797        2.608       16.733       -7.996       -0.028        0.373        0.058        0.043       -0.143
##        2.524        1.206        0.866       67.066       44.006       56.516       -9.239       -7.846       59.577       -0.319       -0.228       -0.166       -0.731        0.547
##       -1.707        0.000        0.000       -2.822       -0.045       -1.912        0.604        2.261       -5.196       11.109        0.000        0.000        4.439        2.410
##        0.000        1.129        0.000        4.742        1.912        2.926        1.826        1.890       -2.279        0.000       10.293        0.000        6.601        0.182
##        0.000        0.000        4.164       -1.313       -0.577       -1.421        2.914        2.676       -6.913        0.000        0.000       10.766        7.967       -0.315
##       -0.070       -0.313       -0.120       -0.339        0.106       -0.181        0.181        0.383       -0.573        0.186        0.251        0.161        1.614       -0.772
##        0.027        0.131        0.037        0.121       -0.035        0.060       -0.061       -0.142        0.215       -0.058       -0.096       -0.060       -0.545        0.271
## 
## Covariance matrix of X (correlations below diagonal):
## 
##       34.083        1.910       -1.209        3.732       -7.305       -0.085       -0.176       -1.045        6.946       -0.287        0.230       -0.458       -0.261        0.509
##        0.056       33.707       -1.025        1.256       -6.446       -1.619        0.887       -1.084        2.984       -0.081        0.083        0.161        0.374       -0.820
##       -0.037       -0.031       31.681        2.346       -1.134        0.885        1.801        1.076        4.881        0.317        0.155       -0.311       -0.397        0.559
##        0.063        0.021        0.041      103.567       66.290       89.536       -4.793       -2.391       67.559        0.403        0.308        0.347       -1.895        4.081
##       -0.112       -0.100       -0.018        0.584      124.354       65.441       -5.364        1.561       42.098        1.029       -0.010        0.291        1.601        2.112
##       -0.002       -0.030        0.017        0.931        0.621       89.284       -4.765       -1.498       57.357        0.603        0.288        0.355       -1.555        3.399
##       -0.007        0.033        0.070       -0.103       -0.105       -0.110       21.102        3.330      -12.774        0.297        0.972        0.546        4.317       -2.071
##       -0.039       -0.041        0.042       -0.051        0.030       -0.034        0.158       21.134       -9.977        0.438        1.652        0.602        4.214       -2.245
##        0.140        0.061        0.102        0.783        0.446        0.716       -0.328       -0.256       71.792       -1.080       -2.388       -1.975       -5.922        3.748
##       -0.026       -0.007        0.030        0.021        0.049        0.034        0.034        0.050       -0.067        3.594       -0.044        0.085        2.635       -0.346
##        0.019        0.007        0.013        0.015        0.000        0.015        0.101        0.172       -0.135       -0.011        4.366        0.009        3.354       -0.706
##       -0.041        0.015       -0.029        0.018        0.014        0.020        0.063        0.069       -0.123        0.024        0.002        3.601        3.218       -0.809
##       -0.010        0.015       -0.016       -0.042        0.032       -0.037        0.212        0.206       -0.157        0.313        0.362        0.382       19.708       -9.614
##        0.027       -0.043        0.030        0.122        0.058        0.109       -0.137       -0.149        0.135       -0.056       -0.103       -0.130       -0.659       10.796

Interpretación:

  • La formación de amistades es más probable con personas del mismo género y el mismo hábito de fumar.

  • Los efectos de tabaquismo de las personas y del ego no son significativos.

  • El estado de fumador aumenta con el tiempo de forma lineal. No se cataloga como efecto significativo en el modelo, porque el poder de detección para este efecto crece en el orden \(O(n)\), siendo \(n\) el número de nodos de la red, y dado que son pocos actores (además de pocos instantes de tiempo observados), el error estándar de la estimación es muy alto.

Bondad de ajuste

Se evaluará la bondad de ajuste sobre del modelo sobre el grado de entrada, propiedad de las redes simuladas que no se incluyó formalmente como predictor en el modelo ajustado.

# Grados de entrada (red del periodo 4)
table(igraph::degree(fr_w4, mode = "in"))
## 
##  1  2  3  4  5  6  7  8  9 10 
##  2  3  7  6  5  4  6  2  1  1
# Evaluar bondad de ajuste
gofi <- sienaGOF(RSmod1, IndegreeDistribution,
                 levls = 1:10, verbose = FALSE, join = TRUE,
                 varName="fr4wav")

# Grafico de bondad de ajuste sobre grados de entrada (1 a 10)
plot(gofi, ylab = "Estadística", 
     main = "Bondad de ajuste de la distribución del grado de entrada")

El modelo reproduce correctamente el grado de popularidad de los actores.

El valor-p de alrededor de 0.8 hace referencia a la ubicación de la distancia de Mahalanobis de la estadística de interés observada, con respecto a la media de la misma en las redes simuladas. Se espera que el valor-p sea cercano a 1, porque eso nos dirá que la distancia de Mahalanobis con respecto a la media es casi 0, es decir, la estadística observada y la media de la misma, dada por las simulaciones, son muy cercanas.

Simulación contra observación

library(igraph)

el <- RSmod1$sims[[500]][[1]][[1]][[3]]
sb <- RSmod1$sims[[500]][[1]][[2]][[3]]

fr_w4_sim <- graph_from_data_frame(el,directed = TRUE)

V(fr_w4_sim)$smoke <- sb
V(fr_w4_sim)$gender <- V(fr_w4)$gender


colors <- c("darkgreen","lightskyblue2")
coord <- layout_with_kk(fr_w4)

# Visualizacion de observada vs. simulada
layout(matrix(c(1, 2,
                3, 3), nrow = 2, byrow = TRUE),
       heights = c(10, 1))

# Gráfico 1
par(mar = c(1,1,3,1))
plot(fr_w4,
     vertex.color = colors[V(fr_w4)$smoke+1],
     vertex.shape = shapes[V(fr_w4)$gender],
     vertex.size = 10,
     main = "Observada - Oleada 4",
     vertex.label = NA,
     edge.arrow.size = 0.2,
     layout = coord)

# Gráfico 2
par(mar = c(1,1,3,1))
plot(fr_w4_sim,
     vertex.color = colors[V(fr_w4_sim)$smoke+1],
     vertex.shape = shapes[V(fr_w4_sim)$gender],
     vertex.size = 10,
     main = "Simulada - Oleada 4",
     vertex.label = NA,
     edge.arrow.size = 0.2,
     layout = coord)

# Panel para la leyenda
par(mar = c(0,0,0,0))
plot.new()

legend("center",
       legend = c("No fumador", "Fumador", "Mujer", "Hombre"),
       col = c("darkgreen", "lightskyblue2", "black", "black"),
       pch = c(16, 16, 1, 0), horiz = TRUE, pt.cex = 1.4, bty = "n")

Comparación de red observada y red simulada

A simple vista la red parece provenir del mismo proceso generador dado por el modelo.

Referencias