Fixed a bug where could not be accessed from a model object, and adjusted one test...
[openmx:openmx.git] / demo / TwoLocusLikelihood.R
1 #
2 #   Copyright 2007-2015 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 #      Hermine Maes -- 2014.11.02 piecewise specification
34 # -----------------------------------------------------------------------------
35
36
37 require(OpenMx)
38 # Load Library
39 # -----------------------------------------------------------------------------
40         
41 # Matrices for allele frequencies, p and s
42 matP         <- mxMatrix( type="Full", nrow=1, ncol=1, 
43                           free=TRUE, values=c(.3333), name="P")
44 matS         <- mxMatrix( type="Full", nrow=1, ncol=1, 
45                           free=TRUE, values=c(.3333), name="S")
46 # Matrix of observed data    
47 obsFreq      <- mxMatrix( type="Full", nrow=4, ncol=1, 
48                           values=c(211,104,39,148), name="ObservedFreqs")
49 matQ         <- mxAlgebra( expression=1-P, name="Q")
50 matT         <- mxAlgebra( expression=1-S, name="T")
51 # Algebra for predicted proportions
52 expFreq      <- mxAlgebra( rbind ((P*P+2*P*Q)*T*T, (Q*Q)*(S*S+2*S*T), 
53                           (P*P+2*P*Q)*(S*S+2*S*T), (Q*Q)*(T*T)), name="ExpectedFreqs")
54 # Algebra for -2logLikelihood
55 m2ll         <- mxAlgebra( expression=-(sum(log(ExpectedFreqs) 
56                           * ObservedFreqs)), name="NegativeLogLikelihood")
57 # User-defined objective
58 funAl        <- mxFitFunctionAlgebra("NegativeLogLikelihood") 
59
60 TwoLocusModel <- mxModel("TwoLocus",
61                           matP, matS, matQ, matT, obsFreq, expFreq, m2ll, funAl)
62
63 # Create an MxModel object
64 # -----------------------------------------------------------------------------
65
66 TwoLocusFit<-mxRun(TwoLocusModel, suppressWarnings=TRUE)
67 TwoLocusFit$matrices
68 TwoLocusFit$algebras
69
70
71 estimates <- c(
72     TwoLocusFit$matrices$P$values, 
73     TwoLocusFit$matrices$S$values, 
74     TwoLocusFit$algebras$NegativeLogLikelihood$result)
75 Mx1Estimates<-c(0.2915,0.1543,647.894)
76
77 omxCheckCloseEnough(estimates,Mx1Estimates,.01)
78 # Compare OpenMx estimates to original Mx's estimates
79 # -----------------------------------------------------------------------------