updates
[openmx:openmx.git] / demo / TwoLocusLikelihood.R
1 #
2 #   Copyright 2007-2009 The OpenMx Project
3 #
4 #   Licensed under the Apache License, Version 2.0 (the "License");
5 #   you may not use this file except in compliance with the License.
6 #   You may obtain a copy of the License at
7
8 #        http://www.apache.org/licenses/LICENSE-2.0
9
10 #   Unless required by applicable law or agreed to in writing, software
11 #   distributed under the License is distributed on an "AS IS" BASIS,
12 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 #   See the License for the specific language governing permissions and
14 #   limitations under the License.
15
16 # -----------------------------------------------------------------------
17 # Program: TwoLocusLikelihood.R  
18 #  Author: Michael Neale
19 #    Date: 11 23 2009 
20 #
21 # Two Locus Likelihood model to estimate allele frequencies
22 #     Bernstein data on ABO blood-groups
23 #     c.f. Edwards, AWF (1972)  Likelihood.  Cambridge Univ Press, pp. 39-41
24 #
25 # Revision History
26 #   Hermine Maes -- 02 22 2010 updated & reformatted
27 # -----------------------------------------------------------------------
28
29
30 require(OpenMx)
31
32 TwoLocusModel <- mxModel("TwoLocus",
33 # Matrices for allele frequencies, p and s
34     mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=c(.3333), name="P"),
35     mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=c(.3333), name="S"),
36 # Matrix of observed data
37     mxMatrix( type="Full", nrow=4, ncol=1, values=c(211,104,39,148), name="ObservedFreqs"),
38 # Algebra for predicted proportions
39     mxAlgebra( expression=1-P, name="Q"),
40     mxAlgebra( expression=1-S, name="T"),
41     mxAlgebra(rbind ((P*P+2*P*Q)*T*T, (Q*Q)*(S*S+2*S*T), (P*P+2*P*Q)*(S*S+2*S*T), (Q*Q)*(T*T)), name="ExpectedFreqs"),
42 # Algebra for -logLikelihood
43     mxAlgebra( expression=-(sum(log(ExpectedFreqs) * ObservedFreqs)), name="NegativeLogLikelihood"),
44 # User-defined objective
45     mxAlgebraObjective("NegativeLogLikelihood")
46 )
47
48 TwoLocusFit<-mxRun(TwoLocusModel)
49 TwoLocusFit@matrices
50 TwoLocusFit@algebras
51
52 # Compare OpenMx estimates to original Mx's estimates
53 estimates <- c(
54     TwoLocusFit@matrices$P@values, 
55     TwoLocusFit@matrices$S@values, 
56     TwoLocusFit@algebras$NegativeLogLikelihood@result)
57 Mx1Estimates<-c(0.2915,0.1543,647.894)
58 omxCheckCloseEnough(estimates,Mx1Estimates,.01)
59