# Fitting a Poisson Latent Space Model for networks

Hi! I’m trying to fit the Latent Space Model (LSM) for an undirected weighted network.

This is the original paper, discussing the model for unweighted networks https://sites.stat.washington.edu/raftery/Research/PDF/hoff2002.pdf

The model, basically, is ` y_ij ~ Poisson(alpha - distance_ij)` i.e. each edge weight is modelled as a Poisson of parameter “some intercept minus the distance between the two extremities of the edge”. I use the squared distance because it’s faster.

The model is not identifiable for the positions in the latent space i.e. unvariant under rotation, translation, reflection. What I do (not included here) is do a Procrustes transform after on the positions. That part works fine.

In the generated quantities, I calculate distances between all positions. Those distances, however, should be identifiable.

My problem is: when the data is not too noisy, the distances are identifiable and that’s good. But when the data is noisy, there are problems in identifying even the distances.

Below is the stan model. Below the stan model are two examples with fake data. In the first example, the data is not too noisy. The distances (generated quantities) converge. In the second example, when the data is noisy, the chains do not converge.

In a third and final example, at the end, I use my real data and nothing converges.

``````data {
int<lower=0> N; // n edges
int<lower=0> M; // n nodes
int<lower=0> D; // D=2
array[N] int<lower=0> y; // edge weight
array[N] int<lower=1> from; // node i
array[N] int<lower=1> to; // node j
}
parameters {
real alpha; // intercept
matrix[M, D] K; // positions
}
transformed parameters {
array[N] real mu;
for (n in 1:N){
mu[n] = exp(alpha-squared_distance(K[from[n],],K[to[n],]));
}

}
model {
// priors
alpha ~ normal(1,3);
to_vector(K) ~ normal(0,3);
// likelihood
for (n in 1:N){
y[n] ~ poisson(mu[n]);
}

}
generated quantities{

vector[N] pred;
vector[N] log_lik;
vector[N] distances;

for (n in 1:N){

pred[n] = poisson_rng(mu[n]);
log_lik[n] = poisson_lpmf(y[n] | mu[n]);
distances[n] = squared_distance(K[from[n],],K[to[n],]);

}
}
``````

Making fake data.

``````library(tidyverse)
library(rstan)

n <- 20
set.seed(222)
x <- rnorm(n,0,1)
y <- rnorm(n,0,1)
fake_df <- tibble(x=x,y=y,id=1:n)
mymat <- matrix(rep(0,n*n),ncol=n)
rownames(mymat) <- 1:n
colnames(mymat) <- 1:n

diag(mymat) <- NA
mymat[upper.tri(mymat)] <- NA

mymat <- as.data.frame(mymat) %>%
mutate(from=row.names(.)) %>%
gather(to,value,-from) %>%
filter(from!=to) %>% filter(complete.cases(.)) %>%
select(-value)
euclidean <- function(a, b) sqrt(sum((a - b)^2))
mymat\$edist <- NA
for (i in 1:nrow(mymat)){
edist <- euclidean(unlist(fake_df[as.numeric(mymat\$from[i]),1:2]),
unlist(fake_df[as.numeric(mymat\$to[i]) ,1:2]))
mymat\$edist[i] <- edist
mymat\$valuenottoonoisy[i] <- rpois(1,exp(3-edist+rnorm(1,0,0.1)))
mymat\$valuenoisy[i] <- rpois(1,exp(3-edist+rnorm(1,0,1.25)))
}
words <- as.character(1:n)
data <- mymat
data\$from_n <- as.numeric(factor(data\$from,words))
data\$to_n <- as.numeric(factor(data\$to,words))

plot(data\$edist,data\$valuenottoonoisy)
plot(data\$edist,data\$valuenoisy)

stan_data <- list(N=nrow(data),
M=length(unique(c(data\$from,data\$to))),
D=2,
y=data\$valuenottoonoisy,
from=data\$from_n,
to=data\$to_n)
``````

Fitting the data to the not too noisy data.

``````model <- stan_model(file = "model.stan")
fit <- sampling(model,data=stan_data, chains = 4, cores = 4, iter = 2000,
control=list(max_treedepth=15))

as.data.frame(summary(fit)\$summary) %>%
filter(grepl("distance",row.names(.))) %>% pull(Rhat) %>%
summary()
``````

Min. 1st Qu. Median Mean 3rd Qu. Max.
0.9991 0.9997 0.9999 1.0000 1.0003 1.0024

That’s good.

``````stan_data\$y <- data\$valuenoisy
fit <- sampling(model,data=stan_data, chains = 4, cores = 4, iter = 2000,
control=list(max_treedepth=15))

as.data.frame(summary(fit)\$summary) %>%
filter(grepl("distance",row.names(.))) %>% pull(Rhat) %>%
summary()
``````

Min. 1st Qu. Median Mean 3rd Qu. Max.
0.9998 1.0604 1.2471 1.5819 1.7158 8.9999

``````stan_data <- read_rds("https://github.com/justinsavoie/personal_website/raw/main/real_data.rds")
fit <- sampling(model,data=stan_data, chains = 4, cores = 4, iter = 2000,
control=list(max_treedepth=15))

as.data.frame(summary(fit)\$summary) %>%
filter(grepl("distance",row.names(.))) %>% pull(Rhat) %>%
summary()
``````

Min. 1st Qu. Median Mean 3rd Qu. Max.
0.9992 1.0334 1.2037 1.7067 1.8082 17.6371