Testing for a significant difference between ML estimates: Likelihood ratio or Wald test?

Question

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?