[R] Bug or Feature? LogLik.nls and non-central F distribution.

Matthew L. Fidler fidler at math.utah.edu
Wed Jan 15 16:57:26 CET 2003


I have a dataset that I am running non-linear regression on via the
following code:


Hill <- function(E0,Em,C50,g,C){
  #
  # Hill is the hill interaction function.
  #
  # E0 Represents the minimum interaction Effect
  #
  # Em Represents the Maximum Interaction Effect
  #
  # C50 represents the concentration at which 50% of the effect occurs.
  #
  # gamma represents the cooperativity of the interaction curve. 
Somtimes called the hill slope.
  #
  # C is the concentration of drug.
  #
  E0 + (Em-E0)*C^g/(C^g+C50^g)
}

short.alif.data <- read.table("shortAlifentanil.txt",header=TRUE);

short.alif.fit <- nls(
        formula = E ~ Hill(0,1,A50,g,A),
        data = short.alif.data,
        start=list(g=1,A50=0.1)
);

Which gives:

> summary(short.alif.fit);

Formula: E ~ Hill(0, 1, A50, g, A)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
g   5.689027   1.367207   4.161 0.025247 *  
A50 0.093200   0.004438  21.000 0.000236 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Residual standard error: 0.08848 on 3 degrees of freedom

Correlation of Parameter Estimates:
         g
A50 0.1761

> logLik(short.alif.fit);
`log Lik.' 6.307024 (df=1)

My problem with this number is that exp(logLik) > 1.  If the error
structure is independent, identically distributed normal (which I
assumed was the case), then the Likelihood function should only give a
number from 0 to 1...  I am running R 1.6.1.  Perhaps this is due to
some strange assumption of non-independent normals???  Maybe something
else???  I assumed that the Likelihood function should be

-n/2*log(2*pi)-n/2*log(d)-n/2

Where d=MLE (biased) estimator of the variance

My next bug observation is the non-central F distribution:

pf(qf(0.05,1,6-3,lower.tail=FALSE),1,6-3,ncp=1.92,lower.tail=FALSE)

This gives a power of a certain linear model that I have done.  However,
the tail is wrong.  To get the right answer you should type

pf(qf(0.05,1,6-3,lower.tail=FALSE),1,6-3,ncp=1.92,lower.tail=TRUE)

If one knows the bug, it is easy to get around, but... it shouldn't
return the wrong number.


-- 
Matthew L. Fidler
fidler at math.utah.edu

421 Wakra Way, Suite 318
Salt Lake City, UT 84108
(801) 581-7125
-------------------------------------------------------------------------
It looked like something resembling white marble, which was
probably what it was: something resembling white marble.
		-- Douglas Adams, "The Hitchhikers Guide to the Galaxy"




More information about the R-help mailing list