## ## Estimering av allelfrekvenser fra ufullstendige blodtypedata. ## ## Simulert data: pA <- .1 pB <- .1 pO <- 1-pA-pB n <- 500 X <- rmultinom(n=1,size=n,prob=c(pA^2+2*pA*pO,pB^2+2*pB*pO,2*pA*pB,pO^2)) ## likelihood funksjon lnL.ABO <- function(p,X) { pA <- p[1] pB <- p[2] pO <- 1-pA-pB nA <- X[1] nB <- X[2] nAB <- X[3] nOO <- X[4] n <- sum(X) lnL <- nA * log(pA^2 + 2*pA*pO) + # AA og AO nB * log(pB^2 + 2*pB*pO) + # BB og BO nAB * log(2*pA*pB) + # AB nOO * log(pO^2) # OO return(-lnL) } ## sme X <- c(44,27,4,88) # Data fra Crow (1986), s. 24 optim(c(.1,.1),lnL.ABO,X=X,method="BFGS")$par optim(c(.1,.1),lnL.ABO2,X=X,method="BFGS")$par sum(X)*optim(c(.1,.1),lnL.ABO,X=X,method="BFGS")$par ## EM-algoritme tar vektor av antall i hver fenotype som innargument ## og returnerer vektor av parameterverdier. EM.ABO <- function(X) { # Observerte data nA <- X[1] nB <- X[2] nAB <- X[3] nOO <- X[4] n <- sum(X) # Startverdier på parameterne pA <- .5 pB <- .1 pO <- 1-pA-pB repeat { # Tar vare på gamle estimat før vi beregner nye... gamle.estimat <- c(pA,pB,pO) # Expectation-steg, # Lager "de data vi mangler" ved å sette disse # lik deres forventning gitt nåværende parameterestimat... nAA <- nA * pA^2 / (pA^2 + 2*pA*pO) # nA * P(AA|AA eller AO) nAO <- nA * 2*pA*pO / (pA^2 + 2*pA*pO) nBB <- nB * pB^2 / (pB^2 + 2*pB*pO) nBO <- nB * 2*pB*pO / (pB^2 + 2*pB*pO) # Maksimerings-steg # vanlig SMEer basert på virkelige og "fiktive data" pA <- (2*nAA + nAO + nAB)/(2*n) pB <- (2*nBB + nBO + nAB)/(2*n) pO <- (2*nOO + nAO + nBO)/(2*n) # eventuellt pO <- 1 - pA - pB # cat(pA, pB, pO, sum(pA,pB,pO), "\n") # Undersøk om endringene i løpet av en iterasjon er tilstrekkelig liten if (sqrt(sum((c(pA,pB,pO)-gamle.estimat)^2)) <.0000001) break() } return(c(pA=pA,pB=pB,pO=pO)) } ## sme beregnet ved bruk av EM-algoritmen EM.ABO(X) ## "Moment-metode" (s. 24 i Crow 1986) ## ## E(nA+nOO)=nA+nOO gir estimat av pB... ## MM.ABO <- function(X) { nA <- X[1] nB <- X[2] nAB <- X[3] nOO <- X[4] n <- sum(X) pB <- 1-sqrt((nA+nOO)/n) pA <- 1-sqrt((nB+nOO)/n) pO <- 1-pA-pB return(c(pA=pA,pB=pB,pO=pO)) }