Updated copyright to 2013 for R/ demo/ models/passing and src/ folders, and also...
[openmx:openmx.git] / demo / TwoLocusLikelihood.R
1 #
2 #   Copyright 2007-2013 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: 2009.11.23 
20 #
21 # ModelType: Locus Likelihood
22 # DataType: Frequencies
23 # Field: Human Behavior Genetics
24 #
25 # Purpose:
26 #      Two Locus Likelihood model to estimate allele frequencies
27 #      Bernstein data on ABO blood-groups
28 #      c.f. Edwards, AWF (1972)  Likelihood.  Cambridge Univ Press, pp. 39-41
29 #
30 # RevisionHistory:
31 #      Hermine Maes -- 2010.02.22 updated & reformatted
32 #      Ross Gore -- 2011.06.15 added Model, Data & Field
33 # -----------------------------------------------------------------------------
34
35
36 require(OpenMx)
37 # Load Library
38 # -----------------------------------------------------------------------------
39
40 TwoLocusModel <- mxModel("TwoLocus",
41         
42     mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=c(.3333), name="P"),
43     mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=c(.3333), name="S"),
44         # Matrices for allele frequencies, p and s
45         # -------------------------------------
46     mxMatrix( type="Full", nrow=4, ncol=1, values=c(211,104,39,148), name="ObservedFreqs"),
47         # Matrix of observed data
48         # -------------------------------------
49     mxAlgebra( expression=1-P, name="Q"),
50     mxAlgebra( expression=1-S, name="T"),
51     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"),
52         # Algebra for predicted proportions
53         # -------------------------------------
54     mxAlgebra( expression=-(sum(log(ExpectedFreqs) * ObservedFreqs)), name="NegativeLogLikelihood"),
55         # Algebra for -logLikelihood
56         # -------------------------------------
57     mxFitFunctionAlgebra("NegativeLogLikelihood")
58         # User-defined objective
59         # ------------------------------------- 
60 )
61 # Create an MxModel object
62 # -----------------------------------------------------------------------------
63
64 TwoLocusFit<-mxRun(TwoLocusModel, suppressWarnings=TRUE)
65 TwoLocusFit@matrices
66 TwoLocusFit@algebras
67
68
69 estimates <- c(
70     TwoLocusFit@matrices$P@values, 
71     TwoLocusFit@matrices$S@values, 
72     TwoLocusFit@algebras$NegativeLogLikelihood@result)
73 Mx1Estimates<-c(0.2915,0.1543,647.894)
74
75 omxCheckCloseEnough(estimates,Mx1Estimates,.01)
76 # Compare OpenMx estimates to original Mx's estimates
77 # -----------------------------------------------------------------------------