www.phidot.org

Hi,

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

I think you want to accomplish this by specifying a design argument to var.components or var.components.reml to take averages. I didn't look to see what your problem was to cause error.

Thank you, var.components and var.components.reml gave much more sensible results.

My motivation for variance components with the whole dataset, then the delta method for a subset of the S-tilde was because of the advice to use at least 10-15 estimates for variance components analysis. The RMark dipper example uses only 6 and I have only 3 to 5, so is there a reason why var.components and var.components.reml are still appropriate?

I still think the delta method on the S-tilde should have worked ok, though obviously it didn't. Is there a reason why this approach should not be used?

Lastly, when var.components and var.components.reml yield similar but different results, should I prefer that from var.components.reml?

Not sure i follow your logic. The stilde estimates were derived from var.components weren't they? So using those with delta method doesn't avoid concern about number of times being used to estimate process variance. My opinion that real is better but others may have thoughts.

Jeff

Ok thank you. I guess having so few estimates is just a limitation in the data.

I was also thinking about finding a robust estimate of the median, rather than the mean.
Is it possible to do this with var.components or var.components.reml?

No there is not.