File: oneLocusLikelihood.R

package info (click to toggle)
r-cran-openmx 2.21.1%2Bdfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: bookworm
  • size: 14,412 kB
  • sloc: cpp: 36,577; ansic: 13,811; fortran: 2,001; sh: 1,440; python: 350; perl: 21; makefile: 5
file content (59 lines) | stat: -rw-r--r-- 2,528 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#
#   Copyright 2007-2020 by the individuals mentioned in the source code history
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#        http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

library(OpenMx)
library(testthat)

if (mxOption(key='Default optimizer') == "CSOLNP") stop("SKIP")

#Single locus Likelihood example
#Author: Mike Neale
#Date: Nov 23 2009

#     Bernstein data on ABO blood-groups
#     c.f. Edwards, AWF (1972)  Likelihood.  Cambridge Univ Press, pp. 39-41
#
onelocus<-mxModel("onelocus", mxFitFunctionAlgebra("NegativeLogLikelihood"),

				mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=c(.3333), name="P"),  # P, freq of allele 1
				mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=c(.3333), name="Q"),  # Q, freq of allele 2
				mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=c(.3333), name="R"),  # R, freq of allele 3
				mxMatrix("Unit", nrow=1, ncol=1, name="I"),                              # 1.0 constant for equality constraint
				mxMatrix("Full", nrow=4, ncol=1, values=c(212,103,39,148),name="ObservedFreqs"), # Data

				mxAlgebra(rbind(P*(P+2*R), Q*(Q+2*R), 2*P*Q, R*R), name="ExpectedFreqs"), # Predicted proportions
				mxAlgebra(-(sum( log(ExpectedFreqs) * ObservedFreqs )), name = "NegativeLogLikelihood"),  # here is -log Likelihood

				mxConstraint(P + Q + R == 1, "equalityConstraint")   # This is the equality constraint
			)

#run the model
run<-mxRun(onelocus)
#run$matrices
#run$algebras

omxCheckCloseEnough(run$algebras$NegativeLogLikelihood$result, 627.028, .1)

# Compare OpenMx estimates to original Mx's estimates
estimates<-c(run$matrices$P$values,run$matrices$Q$values,run$matrices$R$values)
Mx1Estimates<-c(0.2945,0.1540,0.5515)
omxCheckCloseEnough(estimates,Mx1Estimates,.01)

onelocus <- mxModel(onelocus,
                    mxConstraint(P + Q + R == 1, "redundant"))
run <- mxRun(onelocus)
omxCheckCloseEnough(run$algebras$NegativeLogLikelihood$result, 627.028, .1)
estimates<-c(run$matrices$P$values,run$matrices$Q$values,run$matrices$R$values)
omxCheckCloseEnough(estimates,Mx1Estimates,.01)