A Pearson residual is defined as: $r_{i}(\theta)=\frac{y_{i}-E(y_{i}|\theta)}{\sqrt{Var(y_{i}|\theta)}}\tag{1}$ Sum of squared standard residuals $X^{2}$ is: $X^{2}=\sum_ir_i^2\tag{2}$ where $y_{i}$ is $i$th random variable ($i$ = 1, 2, ..., N). $\theta$ is also a random variable since it's about Bayesian inference, but this particular question is about "given $\theta$", so I think we can think $\theta$ not "random". N is the number of dependent variables y.

In a binomial model y[i] ~ dbin(theta, n[i]), it is said that

for known $\theta$, we might expect $X^{2}$ to be around N.

For an example in Appendix A, we can expect $X^{2}$ to be 12, given the known posterior $\theta$ = 0.12 (median value). The prior distribution of $\theta$ is uniform. $i$ has the value of 1, 2, 3, ..., 12 because there are twelve hospitals.

The problem is:

I can prove $X^{2}\approx 12$ through simulation in R, as shown in Appendix B. However, I do not know how to prove mathmatically:

$E(X^{2})=N$

Appendix A:

The Table below shows data representing the mortality rates from 12 English hospitals carrarying out heart surgery on children between 1991 and 1995.

Hospital Operations $n_{i}$ Deaths $y_{i}$

Bristol  143 41
Leicester 187 25 
Leeds 323 24
Oxford 122 23
Guys 164 25
Liverpool 405 42
Southampton 239 24
Great Ormond St 482 53
Newcastle 195 26
Harefield 177 25 
Birmingham 581 58
Brompton 301 31

We can use WinBUGS to calculate the posterior distribution of theta, whose median is 0.12.

# WinBUGS code
model {
    for (i in 1:12){
        y[i] ~ dbin(theta, n[i])
    }
    theta ~ dunif(0, 1)
}

# Data:
list(y=c(41,25,24,23,25,42,24,53,26,25,58,31),
     n=c(143,187,323,122,164,405,239,482,195,177,581,301))

Appendix B:

R codes are used to prove that in the example in Appendix A, for known $\theta$ (i.e., 0.12) we might expect sum of squared standard residuals to be total number of hospitals (i.e., 12).

# number of hospital
N_hosp <- 12  
# theta's value is from posterior mode to 
# check what "we might expect" for known theta
theta <- 0.12 
# number of operations in ith hospital
n=c(143,187,323,122,164,405,239,482,195,177,581,301)
# number of observed deaths in ith hospital
y=c(41,25,24,23,25,42,24,53,26,25,58,31)
# expected y for known theta
y_expect <- array(dim = c(N_hosp))
# Pearson residual
res <- array(dim = c(N_hosp))
res2 <- array(dim = c(N_hosp))

# sample from Binomial to simulate deaths that "we might expect"
# calculate sum of squared stand. residuals
cal_SRes2 <- function(){
  for (i in 1:N_hosp) {
  y_expect[i] <- rbinom(1,n[i],theta)
  res[i] <- (y_expect[i] - n[i]*theta)/sqrt(n[i]*theta*(1-theta))
  res2[i] <- res[i]*res[i]
}
# sum of squared residuals
  sum(res2)
}

# run 50000 times to see if "sum(res2)" is around 12
tmp <- 0
for (j in 1:50000) {
  tmp <- tmp + cal_SRes2()
}
avg_sum_resids <- (1/50000)*tmp  # it is 12, proved