Homecooked solution? [🇷 for BE/BA]

posted by ElMaestro  – Denmark, 2017-08-30 01:18 (2865 d 10:30 ago) – Posting: # 17753
Views: 8,887

Hi all,

since no standard function exists which helps me along, I wrote the stuff below.
The structure of the input file is like in the appendices for Schütz, Labes, Fuglsang's paper in the AAPS Journal.
I am sure the code police will tell me this isn't efficient :-D because the same could be achieved with fewer code lines or some apply family calls. I am sure that is correct. Does the code achieve its goal properly, though? That is a much more interesting question. Not at all throughly tested and certainly not validated.
The purpose is to allow some CI info where we don't ask for "T/R" but perhaps for "R/T" or for "B/D" or "Foo/Bar", whatever.




Extract.Hot.Info=function(TreatmentRatio, Z)
{
  i=regexpr("/", TreatmentRatio)
  T1=substring(TreatmentRatio, 1, i-1)
  T2=substring(TreatmentRatio, i+1, nchar(TreatmentRatio)) ##heck!
  for (i in 1:nrow(Z))
   {
    if (as.character(Z[i,1]) == T1) E1=as.numeric(Z[i,2])
    if (as.character(Z[i,1]) == T2) E2=as.numeric(Z[i,2])
   }
  return(E1-E2)
}

Get.CI.BE=function(FName,  alpha, TreatmentRatio, Stage.Term=F)
{
  A=read.table(FName, header=T) ##Note Var is untransformed
  if (Stage.Term)
   M=lm(log(Var)~factor(Trt)+factor(Seq)+factor(Per)+factor(Subj)+factor(Stg), data=A)
  else
   M=lm(log(Var)~factor(Trt)+factor(Seq)+factor(Per)+factor(Subj), data=A)
 
  df=M$df.residual
  Sc=summary(M)$coefficients
  Trt.Lvls=levels(A$Trt)
  Eff.Lvls=rep(0, length(Trt.Lvls)) ##pseudo effect levels. Usefyl for for extr. of differences
 
  for (i in 1:nrow(Sc))
  for (j in 1:length(Trt.Lvls))
   {
    x=paste("factor(Trt)", as.character(Trt.Lvls[j]), sep="")
    if (rownames(Sc)[i]==x) { Eff.Lvls[j]=Sc[i,1]; SEdiffs=Sc[i,2]}
   }
  Z=data.frame(Trt.Lvls, Eff.Lvls)
  Diff.LogScale=Extract.Hot.Info(TreatmentRatio,Z)
  logL=Diff.LogScale-SEdiffs*qt(p=1-alpha, df=df)
  logU=Diff.LogScale+SEdiffs*qt(p=1-alpha, df=df)
  L=list(PE=exp(Diff.LogScale), Lower=exp(logL), Upper=exp(logU))
  return(L)
}
## call it like Get.CI.BE("Cojones.txt",  0.05, "T/R", F)


Pass or fail!
ElMaestro

Complete thread:

UA Flag
Activity
 Admin contact
23,427 posts in 4,929 threads, 1,678 registered users;
78 visitors (0 registered, 78 guests [including 44 identified bots]).
Forum time: 11:48 CEST (Europe/Vienna)

Half the harm that is done in this world
Is due to people who want to feel important.    T. S. Eliot

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