Bioequivalence and Bioavailability Forum • R code

R code [Regulatives / Guidelines]

posted by ElMaestro – Denmark, 2010-02-25 15:22 (5967 d 09:23 ago) – Posting: # 4817
Views: 6,762

Hmmm, this works on yours truly's machine:


Lund <- function (Linmo, alpha, Verbose)

{

  ## Linmo is the linear model object

  ## alpha is self explanatory

  ## Verbose: set to 1 for more info, otherwise 0



  q=qr(model.matrix(Linmo))$rank

  n=length(residuals(Linmo))

  FLAG=0

  t <- qt((1-alpha/n/2),n-q-1)           # code from d_labes

  isr2 <- (n-q)*(t^2)/(t^2+(n-q)-1)      # code from d_labes

  lcv=sqrt(isr2)                         # code from d_labes

  for (i in 1:n)

   if (abs(rstandard(Linmo)[i])>lcv)

   {

    FLAG=1

    if (Verbose==1) cat ("Stud. residual # ", i, " is ", rstandard(Linmo)[i], ", lcv is ", lcv, "\n")

   }

  if (FLAG==1) cat ("This datafit is not without outliers.\n")

}

EM.

Example - taken from Lund's paper.
y=c(64,60,71,61,54,77,81,93,93,51,76,96,77,93,95,54,168,99)
x1=c(0.4, 0.4, 3.1, 0.6, 4.7, 1.7, 9.4, 10.1, 11.6, 12.6, 10.9, 23.1, 23.1, 21.6, 23.1, 1.9, 26.8, 29.9)
x2=c(53,23,19,34,24,65,44,31,29,58,37,46,50,44,56,36,58,51)

Lm2=lm(y~x1+x2)
Lund(Lm2, 0.01,1)

Complete thread:

Canada: Draft Guidances published Helmut 2010-01-28 00:52
- Outlier earlybird 2010-02-09 08:35
  - Outlier (tests?) Helmut 2010-02-09 14:55
    - Outlier (tests?) earlybird 2010-02-09 16:45
      - Outlier (tests?) Helmut 2010-02-09 17:05
        
        Outlier (tests?) earlybird 2010-02-09 17:18
        Outlier "im Quadrat" d_labes 2010-02-10 13:51
        
        Hotelling's T² Helmut 2010-02-11 17:56
    - Lund! d_labes 2010-02-10 13:36
      - Lund! Helmut 2010-02-12 20:04
        
        intercept ElMaestro 2010-02-12 20:29
        
        intercept Helmut 2010-02-13 00:03
        
        intercept ElMaestro 2010-02-13 01:04
        
        Lund's critical values, technically d_labes 2010-02-22 12:52
        
        R codeElMaestro 2010-02-25 14:22
        
        R code quest d_labes 2010-02-25 16:34
        
        R code quest ElMaestro 2010-02-25 17:06