Delta Method et Var/Cov Matrix

[B.K. Sandercock]

Okay, dumb error on my part for the last example. Revised code as follows. If I take the reported values for mass2 instead, I get the correct estimate for phihat but the variance is still slightly different. For mass = 110, the estimated variance is 9.488912e-05 vs. 7.387e-05 reported in the manual. The small difference might be a rounding error or perhaps a lack of fit to the first-order Taylor series?

# Page B-20 Example 4 variance of back-transformed estimates - a bit harder still
# phihat = exp(-1*(int+beta1+beta2))/(1+exp(-1*(int+beta1+beta2)))
mass = 110; mass.mean = 109.97; mass.SD = 24.79
mass2 = mass^2; mass2.mean = 12707.46; mass2.SD = 5532.03
mass.z = (mass - mass.mean)/mass.SD; mass2.z = (mass2 - mass2.mean)/mass2.SD
par.name = c("int", "beta1", "beta2")
par.mean = c(int= 0.256732, beta1=1.1750358*mass.z, beta2=-1.0554864*mass2.z)
par.varcovar = matrix(c(0.0009006921, -0.0004109710, 0.0003662359,
-0.0004109710, 0.0373887267, -0.0364250288,
0.0003662359, -0.0364250288, 0.0362776933), nrow = 3, byrow = T)
(phihat=plogis(as.numeric(sum(par.mean))))
(phihat.var = deltavar(exp(-1*(int+beta1+beta2))/(1+exp(-1*(int+beta1+beta2))), meanval=par.mean, vars=par.name, Sigma=par.varcovar, verbose=T))
(phihat.SE = sqrt(phihat.var))
# check answers, mass=110, phihat=0.5925, var(phihat) = 0.00007387, SE(phihat) = 0.00860
# check answers, mass=120, phihat=(not given), var(phihat) = 0.000074214, SE(phihat) = 0.008615

[cooch]

par.mean = c(int= 0.256732, beta1=1.1750358*mass.z, beta2=-1.0554864*mass2.z)

As per p. 20 in the Appendix, should be

par.mean = c(int= 0.2567333, beta1=1.1750545*mass.z, beta2=-1.0555046*mass2.z)

Small differences, but in this game, they all add up. Haven't tested...

par.varcovar = matrix(c(0.0009006921, -0.0004109710, 0.0003662359,
-0.0004109710, 0.0373887267, -0.0364250288,
0.0003662359, -0.0364250288, 0.0362776933), nrow = 3, byrow = T)

Here is the VC matrix output from MARK - again, close, but not identical. What you have is what is in 'the book'.

0.0009006961 -0.0004109935 0.0003662581
-0.0004109935 0.0373907263 -0.0364270112
0.0003662581 -0.0364270112 0.0362796588


I'll check more when I have a chance.

[B.K. Sandercock]

Can I trouble you to double-check that you have posted the most current version of Appendix B? The numbers that I used for the coefficients and the var-covar matrix were taken from pages B-20 and B-22, and are the same in the individual chapter, full book, and rotated version of the book. The revised numbers that you posted above presumably explain the difference in estimates of var(phi), but do not match the text of the Appendix for any of the PDFs that I downloaded today. Thanks.

[cooch]

Okay, dumb error on my part for the last example. Revised code as follows. If I take the reported values for mass2 instead, I get the correct estimate for phihat but the variance is still slightly different. For mass = 110, the estimated variance is 9.488912e-05 vs. 7.387e-05 reported in the manual. The small difference might be a rounding error or perhaps a lack of fit to the first-order Taylor series?

I double-checked the calculations, in 2 different CAS systems (Maple, Maxima). They all gave *exactly* the same answer (out to 10+ decimal places), which is in the version the appendix I just posted a few minutes ago (not substantially different from the version Brett was looking at).

The correct estimate of the variance based on a first order approximation is definitely 0.00073866, giving an estimated SE of 0.0085945, which is *exactly* what is reported by MARK.

So, either

1\ Bolker's code has some problems

2\ Brett's code has some problems

3\ Bolker's code is using a higher-order approximation.

I'd say of the 3, possibility (3) seems pretty unlikely, but I can't be sure (seems unlikely since Brett reports that it matches up fine with the other examples in the Appendix, all of which were generated using the first-order approximation). A true test would be to try it against an example where we 'know' the result for a second-order fit. See B.3.1 for an example of same. (I'm afraid I don't have time to test it). If Brett's code is correct (ibid, note that the numbers if Appendix B have been updated to match what MARK now reports for this data set) then that leaves Bolker's code as the prime suspect. To some degree, this is why I prefer a CAS approach - may take longer (although I have it pretty well automated at my own end), but I can see *exactly* what is happening at each step.

[cooch]

Can I trouble you to double-check that you have posted the most current version of Appendix B? The numbers that I used for the coefficients and the var-covar matrix were taken from pages B-20 and B-22, and are the same in the individual chapter, full book, and rotated version of the book. The revised numbers that you posted above presumably explain the difference in estimates of var(phi), but do not match the text of the Appendix for any of the PDFs that I downloaded today. Thanks.

Just did -- see my last post, which crossed yours over the net by about 1 minute.

[egc]

Because this thread has more to do with the 'Delta method generally, and only superficially to E-SURGE (or anything particularly 'software-specific'), I'm going to move this entire thread over to 'General Questions'.

[B.K. Sandercock]

I downloaded the updated Appendix B - thanks for posting. Anybody else doing this, note that you have to click the link for the PDF and then click Reload to get the new version, otherwise you might just retrieve the old version out of your cache. The var-covar matrix on B-22 has the updated numbers, but the beta coefficients on B-20 are still the old numbers.

Here is an R script with the revised numbers.
# Revised numbers
# Page B-20 Example 4 variance of back-transformed estimates - a bit harder still
# phihat = exp(-1*(int+beta1+beta2))/(1+exp(-1*(int+beta1+beta2)))
mass = 110; mass.mean = 109.97; mass.SD = 24.79
mass2 = mass^2; mass2.mean = 12707.46; mass2.SD = 5532.03
mass.z = (mass - mass.mean)/mass.SD; mass2.z = (mass2 - mass2.mean)/mass2.SD
par.name = c("int", "beta1", "beta2")
par.mean = c(int= 0.2567333, beta1=1.1750545*mass.z, beta2=-1.0555046*mass2.z)
par.varcovar = matrix(c(0.0009006961, -0.0004109935, 0.0003662581,
-0.0004109935, 0.0373907263, -0.0364270112,
0.0003662581, -0.0364270112, 0.0362796588), nrow = 3, byrow = T)
(phihat=plogis(as.numeric(sum(par.mean))))
(phihat.var = deltavar(exp(-1*(int+beta1+beta2))/(1+exp(-1*(int+beta1+beta2))), meanval=par.mean, vars=par.name, Sigma=par.varcovar))
(phihat.SE = sqrt(phihat.var))
# check answers, mass=110, phihat=0.5925, var(phihat) = 0.00007387, SE(phihat) = 0.00860
# check answers, mass=120, phihat=(not given), var(phihat) = 0.000074214, SE(phihat) = 0.008615

Unfortunately, the function still gives slightly different answers versus Appendix B.
phihat
[1] 0.5924384
phihat.var
[1] 9.488934e-05
phihat.SE
[1] 0.009741116

I'll follow Evan's advice and try the calculations in a different platform like Matlab. Thanks for the help with trouble-shooting.

Brett.

[cooch]

I downloaded the updated Appendix B - thanks for posting. Anybody else doing this, note that you have to click the link for the PDF and then click Reload to get the new version, otherwise you might just retrieve the old version out of your cache.

You should have an add-on for Firefox or Chrome that does this automatically anyway (i.e., clears your cache when you turn off your browser). Not only does this save you the hassle of not knowing if you have the most current version of something, but it is a 'best practice' for basic security.

The var-covar matrix on B-22 has the updated numbers, but the beta coefficients on B-20 are still the old numbers.

Yeah, but close enough. I had rather assumed people would not use the values written in the book - which are rounded to make the formatting reasonable - but would actually run the analysis, and cut-and-paste everything from their output, ensuring full precision. I see what Brett is trying to do - kindly provide a script with all the numbers so that people using it don't have to do any work at all. Well, given that Brett is being so kind, I suggest he take it one step further, run the analysis himself, and copy and paste the estimates (betas, VC) into his script.

Moreover, be advised that the numbers you get might differ >5th decimal place, depending on the optimization algorithm you use, and whether you're running 32- or 64-bit MARK. I occasionally flip back and forth between the default procedures, and simulated annealing. Depending on the problem, you might notice the numbers differ somewhat, especially for the VC matrix. Generally these differences aren't large enough to be important to the final result.

Here is an R script with the revised numbers.
# Revised numbers
# Page B-20 Example 4 variance of back-transformed estimates - a bit harder still
# phihat = exp(-1*(int+beta1+beta2))/(1+exp(-1*(int+beta1+beta2)))
mass = 110; mass.mean = 109.97; mass.SD = 24.79
mass2 = mass^2; mass2.mean = 12707.46; mass2.SD = 5532.03
mass.z = (mass - mass.mean)/mass.SD; mass2.z = (mass2 - mass2.mean)/mass2.SD
par.name = c("int", "beta1", "beta2")
par.mean = c(int= 0.2567333, beta1=1.1750545*mass.z, beta2=-1.0555046*mass2.z)
par.varcovar = matrix(c(0.0009006961, -0.0004109935, 0.0003662581,
-0.0004109935, 0.0373907263, -0.0364270112,
0.0003662581, -0.0364270112, 0.0362796588), nrow = 3, byrow = T)
(phihat=plogis(as.numeric(sum(par.mean))))
(phihat.var = deltavar(exp(-1*(int+beta1+beta2))/(1+exp(-1*(int+beta1+beta2))), meanval=par.mean, vars=par.name, Sigma=par.varcovar))
(phihat.SE = sqrt(phihat.var))
# check answers, mass=110, phihat=0.5925, var(phihat) = 0.00007387, SE(phihat) = 0.00860
# check answers, mass=120, phihat=(not given), var(phihat) = 0.000074214, SE(phihat) = 0.008615


Unfortunately, the function still gives slightly different answers versus Appendix B.
phihat
[1] 0.5924384
phihat.var
[1] 9.488934e-05
phihat.SE
[1] 0.009741116

I'll follow Evan's advice and try the calculations in a different platform like Matlab. Thanks for the help with trouble-shooting.


Brett.

Thats a waste of time, and implies that I have it wrong somehow. I don't. Two different CAS give *exactly* the result reported by MARK. I could run it in MATLAB, or Octave, or Derive, or SciPy, but that would be redundant verging on silly. I would bet significant sums of $$$ they'd give 'exactly' the same answer.

Basically, it comes down to this -- if you have all the right values (betas, VC matrix), and are sure your script is right, then Bolker's function has a problem. Suggesting otherwise implies that not only is MARK wrong, but very high-end CAS environments.

[B.K. Sandercock]

I'm convinced there is problem with the deltavar() function too. I had not appreciated that the var-covar matrix would differ slightly among different runs with different optimization routines, so trying to match up numbers beyond the 5th decimal place has little effect on the calculations. My goal here was to take the R function for a test run and my questions were aimed at trying to reconcile the difference in results. No slight intended by suggesting Matlab, just the only CAS software I've used previously. Had not heard of Maxima before, will have to read more, thanks for the tip.

[cooch]

I'm convinced there is problem with the deltavar() function too. I had not appreciated that the var-covar matrix would differ slightly among different runs with different optimization routines, so trying to match up numbers beyond the 5th decimal place has little effect on the calculations. My goal here was to take the R function for a test run and my questions were aimed at trying to reconcile the difference in results. No slight intended by suggesting Matlab, just the only CAS software I've used previously. Had not heard of Maxima before, will have to read more, thanks for the tip.

No worries. Before you find fault with Bolker's code, try it with full-precision parameters, VC, and make sure your Jacobian vectors match with those in the book (more or less).

Maxima is a free CAS - it isn't perfect, but (like R) it shares the attribute of being free. As such, I use it a lot in classes, since I can ask (require) students to use it (since it doesn't cost them anything). It has a few quirks, but so do the big commercial packages (Maple, Mathematica). Maxima has also been ported to every platform I can think of, so that isn't a limit.

In some ways, MATLAB is the best overall option since it licenses the symbolic toolbox from Maple (meaning, you can do all the linear algebera, plus all the 'other math') in one environment. And the add-on toolboxes for MATLAB are first rate. But, I stopped using it when I realized I could use R for 95% of my 'linear algebra work', and Maxima for everything else. And full-blown MATLAB is wicked expensive (although our local site license is pretty reasonable). I keep *all* of the different software applications around to serve as checks on each other (which is why I'm pretty sure that the Delta results for the problem in question in this thread are correct).

At some point tonite, I'll go through and make sure all pages for that example have the same parameter estimates in place.

Finally, as an example of the difference optimization method can make...consider the following complete VC matrices (for all 4 parameters: phi and p)

Using the default optimization, default starting values:

Code: Select all

0.0009006977   -0.0004110021   0.0003662665   -0.0006544039   
-0.0004110021   0.0373919352   -0.0364281943   0.0000612919   
0.0003662665   -0.0364281943   0.0362808169   -0.0000567379   
-0.0006544039   0.0000612919   -0.0000567379   0.0029274838   


Using simulated annealing, default starting values

Code: Select all

0.0009006949   -0.0004109807   0.0003662452   -0.0006544000   
-0.0004109807   0.0373889730   -0.0364252631   0.0000612911   
0.0003662452   -0.0364252631   0.0362779160   -0.0000567368   
-0.0006544000   0.0000612911   -0.0000567368   0.0029274768   


Using simulated annealing, default method estimates as starting values

Code: Select all

0.0009006966   -0.0004110193   0.0003662835   -0.0006544013   
-0.0004110193   0.0373926588   -0.0364289097   0.0000612978   
0.0003662835   -0.0364289097   0.0362815245   -0.0000567437   
-0.0006544013   0.0000612978   -0.0000567437   0.0029274687   


All 3 are slightly different, albeit fairly 'far out to the right', if you catch my meaning.

However, if you are 'consistent' (i.e., use beta and VC estimates for a given 'method'), then the value you calculate for Var/SE will match the one MARK gives you for that 'method'.