greg
New Member
Posts: 1
|
Post by greg on Apr 28, 2014 0:16:26 GMT
I have been working on extending the gcta approach but am getting slightly different results from the gcta program on the same data with a continuous phenotype.
The easiest way to see the difference is to take the special case where there are no fixed effects, the phenotype is standardized, and the model is V(A) = 0. The log likelihood reduces to -.5*(n*log(Ve) + t(y)*y/Ve). Because the phenotype is standardized, t(y)*y should be close to n and Ve close to 1, so we expect that log(L) will be close to -.5n. Here is that function for my test data: > n <- length(y) > n [1] 1646 > tyty <- t(y) %*% y > tyty [,1] [1,] 1645 > Ve <- seq(from=.99, to=1.01, by=.005) > fVe <- function(Ve) -.5*(n*log(Ve) + tyty/Ve) > fVe(Ve) [1] -822.5367 -822.5078 -822.5000 -822.5127 -822.5456
gcta, however, gives log(L) = -826.203 for this model.
It is not the case that gcta is adding a constant to log(L) because I get a different LRT when I compare this model with the one in which V(A) is free. Also, pruning cases is not an issue because all off diagonal values of the GRM are less than .025.
Puzzled.
|
|