0

I'm having trouble estimating the standard errors from the predicted probabilities from a ERGM model, to calculate a confidence interval. Getting the predicted probabilities is not a problem, but I want to get a sense of the uncertainty surrounding the predictions.

Below is a reproduceable example based on the data set of marriage and business ties among Renaissance Florentine families.

library(statnet)
data(flo)
flomarriage <- network(flo,directed=FALSE)
flomarriage

flomarriage %v% "wealth" <- c(10,36,27,146,55,44,20,8,42,103,48,49,10,48,32,3)
flomarriage
gest <- ergm(flomarriage ~ edges + 
               absdiff("wealth"))
summary(gest)

plogis(coef(gest)[['edges']] + coef(gest)[['absdiff.wealth']]*10)

Based on the model, it is estimated that a wealth difference of 10 corresponds to a 0.182511 probability of a tie. My first question is, is this a correct interpretation? And my second question is, how does one calculate the standard error of this probability?

Christophe Lesschaeve
  • 13
  • 2

2 Answers2

1

This is the correct interpretation for a dyad-independent model such as this one. For a dyad-dependent model, it would be the conditional probability given the rest of the network.

You can obtain the standard error of the prediction on the logit scale by rewriting your last line as a dot product of a weight vector and the coefficient vector:

eta <- sum(c(1,10)*coef(gest))
plogis(eta)

then, since vcov(gest) is the covariance matrix of parameter estimates, and using the variance formulas (e.g., http://www.stat.cmu.edu/~larry/=stat401/lecture-13.pdf),

(var.eta <- c(1,10)%*%vcov(gest)%*%c(1,10))

You can then get the variance (and the standard error) of the predicted probability using the Delta Method (e.g., https://blog.methodsconsultants.com/posts/delta-method-standard-errors/). However, unless you require the standard error as such, as opposed to a confidence interval, my advice would be to first calculate the interval for eta (i.e., eta + qnorm(c(0.025,0.975))*sqrt(var.eta)), then call plogis() on the endpoints.

Pavel N. Krivitsky
  • 66
  • 3
1

Really very helpful answer! Many thanks for it. Here's just a small, almost cosmetic note about the suggested formula: 1) the upper bound in qnorm should be 0.975, not 0.0975, and 2) use matrix multiplication before sqrt --> eta + qnorm(c(0.025,0.975))%*%sqrt(var.eta).

Milchfilm
  • 11
  • 1