"By popular demand": likelihood [R for BE/BA]

posted by ElMaestro  – Denmark, 2020-07-24 10:07 (187 d 05:40 ago) – Posting: # 21783
Views: 8,334

Hi all,

to make the objective function return the likelihood, use e.g.:

Obj.F12=function(Pars)
##returns the restricted log likelihood.
##apart from k, this is line 3 of page 10 of:
##http://people.csail.mit.edu/xiuming/docs/tutorials/reml.pdf
##k is invariant with parameter estimates and variance components.
{
  CovM=Create.CovM(Pars)
  k=  -((length(y)-ncol(X))/2)*log(2*pi)
  A= -0.5*log( det(CovM))
  B= -0.5*log( det(t(X) %*% solve(CovM) %*% X))
  est.b = solve(t(X) %*% solve(CovM) %*% X) %*% t(X) %*% solve(CovM) %*% y
  tmp= y - X %*% est.b
  C=-0.5 *(t(tmp) %*% solve(CovM) %*% tmp)
  return(k+A+B+C)
}


The optimised object (F in my original code) will hold the optimised likelihood as F$value.
Opportunity: Will run a tiny bit faster if k is public, and if solve(CovM) is done a single time.


For example at reltol=1e-8 I get with EMA dataset II a value of 15.33728.

This matches exactly what SAS reports, see post # 21670 above.

Pass or fail!
ElMaestro

Complete thread:

Activity
 Admin contact
21,316 posts in 4,446 threads, 1,489 registered users;
online 4 (0 registered, 4 guests [including 2 identified bots]).
Forum time: Wednesday 14:47 CET (Europe/Vienna)

Nothing fails like success because you do not learn anything from it.
The only thing we ever learn from is failure.
Success only confirms our superstitions.    Kenneth E. Boulding

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5