11
$\begingroup$

I am trying to test whether or not there is a significant difference between maximum likelihood estimates of two genetic parameters (selection and dominance) across two environments with genotype data from a cross. There are three genotype classes, C11, C12, and C22, and I am looking to maximize the likelihood below:

llhood2 <- function(theta, c, c11,c12,c22){

s<-theta[1]
d<-theta[2]

P11 = (1 - ((s * (1 - c) * (1 - c + (2 * d * c)))) / (4 - (s * (1 + (2 * d))))
P22 = (1 - (s * c * ((2 * d * (1 - c)) + c))) / (4 - (s * (1 + (2 * d))))
P12 = (1 - P11 - P22)
logL = c11*log(P11) + c12*log(P12) + c22*log(P22)
-logL   
}

C is fixed depending on which genetic marker you are using, and c11, c12, and c22 are the data, representing genotypes AA, AB, and BB respectively. So I am jointly estimating s and d.

Using the optim function in R and plugging in genotype data from each environment, I am able to get a point estimate for s and d in one environment (s1 and d1) and the variance of these estimates from the hessian, then I get parameter estimates using genotype data from the 2nd environment,(s2 and d2, and thier associated variances). So, two separate fittings of the model give me s1 var_s_1 and s2 var_s_2, etc. The Code for optim is below:

fit1 <- optim(inits, llhood, method="L-BFGS-B", hessian=T, lower=lower, upper=upper,
c=c,c11=c11,c12=c12,c22=c22, control=list(trace=1))

fit2 <- optim(inits, llhood, method="L-BFGS-B", hessian=T, lower=lower, upper=upper,
c=c,c11=n11,c12=n12,c22=n22, control=list(trace=1)

I then want to know if s1 is significantly different from s2 (likewise for d1 and d2). I thought this would be a straight forward test, but I have had a hard time finding documentation on how to do this properly, or most elegantly. As a crude approximation, I used a simple two-sample T test, but Im pretty sure that it is not quite right...

Ive also considered a likelihood ratio test where I test whether the difference between the two estimates of s1 and s2, is not different from zero (reduced model) vs different from zero. But Im not sure how to set this up given my likelihood equations.

Any specific help would be much appreciated.

UPDATE: Here is some attempted code for the restricted and unrestricted likelihoods which will be used to perform a likelihood ratio test. First the restricted (null) where s and ss are equal- (c11-c22 are the data for env.1 , n11-n12 for environment 2)

llhoodnull <- function(theta,c,c11,c12,c22,n11,n12,n22){

s<-theta[1] 
d<-theta[2]
ss=s
dd=d


P11 = (1 - (s * (1 - c) * (1 - c + (2 * d * c)))) / (4 - (s * (1 + (2 * d))))
P22 = (1 - (s * c * ((2 * d * (1 - c)) + c))) / (4 - (s * (1 + (2 * d))))
P12 = (1 - P11 - P22)
b11 = (1 - (ss * (1 - c) * (1 - c + (2 * dd * c)))) / (4 - (ss * (1 + (2 * dd))))
b22 = (1 - (ss * c * ((2 * dd * (1 - c)) + c))) / (4 - (ss * (1 + (2 * dd))))
b12 = (1 - b11 - b22)

logL = c11*log(P11) + c12*log(P12) + c22*log(P22)+n11*log(b11) + n12*log(b12) + 
    n22*log(b22)
-logL

}

I think that I may also need to fix the d so that I am only testing for the likelihood of s.

And the unrestricted:

llhoodb <- function(theta,c,c11,c12,c22,n11,n12,n22){

s<-theta[1] 
d<-theta[2]
ss<-theta[3]
dd<-theta[4]


P11 = (1 - (s * (1 - c) * (1 - c + (2 * d * c)))) / (4 - (s * (1 + (2 * d))))
P22 = (1 - (s * c * ((2 * d * (1 - c)) + c))) / (4 - (s * (1 + (2 * d))))
P12 = (1 - P11 - P22)
b11 = (1 - (ss * (1 - c) * (1 - c + (2 * dd * c)))) / (4 - (ss * (1 + (2 * dd))))
b22 = (1 - (ss * c * ((2 * dd * (1 - c)) + c))) / (4 - (ss * (1 + (2 * dd))))
b12 = (1 - b11 - b22)

logL = c11*log(P11) + c12*log(P12) + c22*log(P22)+n11*log(b11) + n12*log(b12) + 
    n22*log(b22)
-logL

}

This one gives me errors when I try to maximize it with optim()...

Any thoughts on these?

$\endgroup$
  • 2
    $\begingroup$ The most straightforward way to do this, I think, is by maximizing the (log) likelihood of both samples in a single call to optim(). Optimize the likelihood when s1=s2, and then see how much better you do without this constraint. Twice the difference in log-likelihood will have a $\chi^2_1$ distribution, for large samples, if the values of s are actually identical. (See any text on likelihood ratio tests) $\endgroup$ – guest Apr 17 '12 at 3:16
  • $\begingroup$ @guest OK. Any suggestion how to code that to a single call to optim()? I have two different sets of genotype data (one for each environment) that Id have to work in. I also have two unknowns, s and d. so I have to do this for s1=s2 and d1=d2. $\endgroup$ – LPP Apr 17 '12 at 14:30
  • $\begingroup$ If I understand @guest correctly, you need two log-likelihood functions: L1(s1, s2, d; data) and L2(s, d; data). L2 is just like the llhood2 function that you already have. And L1, if your log-likelihood is an aggregation of terms, I think you can do this: L1(s1, s2, d; data)=llhood2(s1, d; subset1)+llhood2(s2, d; subset2), so you only need a wrapper function that splits data into subset1 and subset2 and then call llhood2. $\endgroup$ – Ernest A Apr 17 '12 at 15:02
  • $\begingroup$ @Ernest ok. In the setup you just described, is L1 the reduced or constrained model (where s1=s2?) -- I tried to setup a new likelihood where s is the same for both sets of data and i sum the likelihoods Lnull(s1=s2,d,data1, data2) and then he unrestricted L2(s1,s2,d,data1 data2) where I estimate each s seperately and sum the resultant likelihoods. But the unrestricted likelihood, L2, returns "needs finite values of fn". $\endgroup$ – LPP Apr 17 '12 at 15:20
  • $\begingroup$ It means that L2 returned a non-finite value (Inf, Nan...) most likely because P11, P22 or P12 were non-positive. On trick is to replace all values that are equal or below zero with a very small positive number (e.g 0.000001) before calculating the logerithms. $\endgroup$ – Ernest A Apr 17 '12 at 16:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.