I am attempting to calculate the average and SE of four model estimates, subject to process variance alone. However, I am getting an error after the delta method. Could you please let me know (i) why the error occurs, and (ii) if my method is correct. At the end I have reproduced the calculation with total variance estimates without any problem.

My model has 16 Psi estimates, so I use Variance Components (random effects, intercept only) on all 16 Betas to obtain the 16 S-tilde which are subject only to process variance.

In order to calculate the average of the first 4 real Psi (subject only to process variance), I see that they can be calculated from the five logit-scaled S-tilde betas (labelled x1 to x5);

1st real psi = exp(x1+x2)/(1+exp(x1+x2))

2nd real psi = exp(x1+x3)/(1+exp(x1+x3))

3rd real psi = exp(x1+x4)/(1+exp(x1+x4))

4th real psi = exp(x1+x5)/(1+exp(x1+x5))

Then their average is;

- Code: Select all
`x1 = 1.480638; x2 = -0.897507; x3 = -1.033954; x4 = -0.673738; x5 = -0.582794`

Average = (exp(x1+x2)/(1+exp(x1+x2)) + exp(x1+x3)/(1+exp(x1+x3)) + exp(x1+x4)/(1+exp(x1+x4)) + exp(x1+x5)/(1+exp(x1+x5)))/4 #0.6633982

And SE from the delta method is;

- Code: Select all
`vcov <- c(0.00208, -0.50326, -0.58330, -0.60619, -0.54785,`

-0.00221, 0.00922, 0.30596, 0.31018, 0.27697,

-0.00218, 0.00241, 0.00673, 0.35633, 0.32541,

-0.00229, 0.00247, 0.00242, 0.00686, 0.32754,

-0.00211, 0.00224, 0.00225, 0.00229, 0.00712)

sigma <-matrix(vcov,nrow=5,ncol=5,byrow=TRUE)

seAv=deltamethod(~(exp(x1+x2)/(1+exp(x1+x2)) + exp(x1+x3)/(1+exp(x1+x3)) + exp(x1+x4)/(1+exp(x1+x4)) + exp(x1+x5)/(1+exp(x1+x5)))/4, c(x1,x2,x3,x4,x5), sigma, ses=T)

#Warning message: In sqrt(diag(new.covar)) : NaNs produced

CHECK WITH TOTAL VARIANCE ESTIMATES

MLE beta estimates:

- Code: Select all
`x1=1.490899; x2=-0.9146954; x3= -1.049563; x4=-0.6876908; x5=-0.5955833`

- Code: Select all
`Av = (exp(x1+x2)/(1+exp(x1+x2)) + exp(x1+x3)/(1+exp(x1+x3)) + exp(x1+x4)/(1+exp(x1+x4)) + exp(x1+x5)/(1+exp(x1+x5)))/4 # 0.6623542`

- Code: Select all
`vcov <- c(0.002135057, -0.002263782, -0.002240721, -0.002351096, -0.002162478,`

-0.002263782, 0.009404006, 0.002471348, 0.002529644, 0.002297931,

-0.002240721, 0.002471348, 0.006842304, 0.002483221, 0.002307514,

-0.002351096, 0.002529644, 0.002483221, 0.006971759, 0.002345244,

-0.002162478, 0.002297931, 0.002307514, 0.002345244, 0.007232530)

sigma <-matrix(vcov,nrow=5,ncol=5,byrow=TRUE)

seAv=deltamethod(~(exp(x1+x2)/(1+exp(x1+x2)) + exp(x1+x3)/(1+exp(x1+x3)) + exp(x1+x4)/(1+exp(x1+x4)) + exp(x1+x5)/(1+exp(x1+x5)))/4, c(x1,x2,x3,x4,x5),sigma,ses=T) # 0.008148399