search

hanchenphd / GMMAT

Generalized linear Mixed Model Association Tests
Other
36 stars 22 forks source link

Saddlepoint approximation in `.pKuonen()` #22

Closed dwoll closed 3 years ago

dwoll commented 3 years ago

First of all, many thanks for making your work available for others to use and study! I'm trying to study implementations of the saddlepoint approximation from Kuonen (1999), and I'd be glad if you could help me with the following questions.

In function .pKuonen(), it seems to me that the following line does not achieve anything. Maybe a return() is missing here? https://github.com/hanchenphd/GMMAT/blob/518f66409fdaf007e563557d38259ba079b3a494/R/SMMAT.R#L1170

In sub-functions kprime0() https://github.com/hanchenphd/GMMAT/blob/518f66409fdaf007e563557d38259ba079b3a494/R/SMMAT.R#L1181

and kpprime0() https://github.com/hanchenphd/GMMAT/blob/518f66409fdaf007e563557d38259ba079b3a494/R/SMMAT.R#L1186

The code is using a different notation compared to Kuonen (1999) p. 930. Kuonen uses sigma^2 in the equations (non-centrality parameter for a chi^2 distribution) whereas you use df+delta. My math skills are quite limited, so it would help me a lot if you could explain the relationship between your implementation and Kuonen's equations. Many thanks in advance!

hanchenphd commented 3 years ago

Thanks for catching line 1170. Yes, this line of code should be wrapped by return() and we will fix it in the next version.

The kprime0 and kpprime0 functions are the first and second derivatives of the cumulant generating function K(ζ). Kuonen used h_i and \sigma_i^2 to denote the degree of freedom and the non-centrality parameter, and they were df and delta in the code. I think there was a missing ζ in the numerator of the second term in Kuonen's K(ζ), and perhaps that could explain the confusion when you take the derivatives.

Best, Han

dwoll commented 3 years ago

Thank you very much for your quick response and for taking the time to help me better understand your implementation!

In Kuonen's PhD thesis equation 3.2.4 (page 21), ζ is also missing from the numerator in the second term of the cumulant generating function K. The first and second derivatives on page 24 also look somewhat different compared to your implementation. That is indeed what is confusing to me. Is there a manuscript where the formulas you are implementing are written down? Thanks again!

hanchenphd commented 3 years ago

I believe it was a typo in Kuonen's paper. The cumulant generating function is basically the log of the MGF, and the MGF of a weighted sum of independent non-central chi-squares is the product of the MGF of each term. See Wikipedia and StackExchange for the MGF of non-central chi-squares.

Kuonen cited Equation 2.3 from Imhof 1961, which is consistent with the links above.

Best, Han

dwoll commented 3 years ago

That makes sense, thank you for the explanation!