ElMaestro
★★★

Denmark,
2013-04-08 22:23
(4419 d 03:21 ago)

Posting: # 10375
Views: 6,970
 

 EMA dataset 1 trouble [RSABE / ABEL]

Post reply
Hi all,

I am fiddling a bit around with scaled BE and trying to use R for same.
Cannot get to 107.11% - 124.89% for EMA dataset 1 / method A. My PE is correct so the trouble I guess is in somewhere in one of these parts:

1. sqrt(0.5*(1/nSeq1 + 1/nSeq2)
2. df and the critical T-value.

Ad 1:
Looks to me like the lm allows effect estimates for 75 subjects; this means subjects with a missing value are not tossed out. So I guess all is fine with 77 in total; seemingly 38 and 39. But I get the wrong CI anyway.

Ad 2:
For df I tried both the df form the full model and from the Ref-only model (let me add: I think it is a valid argument that the df from the latter should be used but I am sure we can have a good fight over that some other time).

Any ideas?
Here's some code bits:
A=read.table("EMAset1b.csv", header=T)
Model1=lm(log(A$DATA)~0+factor(A$FORMULATION)+
factor(A$PERIOD)+factor(A$SEQUENCE)+factor(A$SUBJECT))
lnPE=(coef(Model1)[2]-coef(Model1)[1]) ## works!
B=subset(A, FORMULATION=='R')
Model2=lm(log(B$DATA)~
factor(B$PERIOD)+factor(B$SEQUENCE)+factor(B$SUBJECT))
s=summary(Model2)$sigma ## works!
CV=sqrt(exp(s*s)-1.0) ## works!
DF=Model1$df.residual ## or use Model2's residual. Discuss some other time.
CritT=qt(1-0.05, DF)

(blah blah code for extraction of n per sequence)

lnU=lnPE+s*CritT*sqrt(0.5*(1/nBABA+1/nABAB)) ##can't get it right.
lnL=lnPE-s*CritT*sqrt(0.5*(1/nBABA+1/nABAB)) ##can't get it right.

Thanks in advance.
—
Pass or fail!
ElMaestro
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2013-04-09 04:37
(4418 d 21:07 ago)

@ ElMaestro
Posting: # 10376
Views: 5,923
 

 complete (RR) only

Post reply
Hi ElMaestro!

❝ I am fiddling a bit around with scaled BE and trying to use R for same.


Welcome to the crippled model!

❝ Ad 1:

❝ Looks to me like the lm allows effect estimates for 75 subjects; this means subjects with a missing value are not tossed out. So I guess all is fine with 77 in total; seemingly 38 and 39.


Wonderful. Phoenix/WinNonlin by default uses all data as well (REML…). Since we have to mimick SAS we have to get rid of subjects with incomplete data (#24, 31, 67, 71). The complete data set contains 73 subjects with two administrations of R.

❝ Ad 2:

❝ For df I tried both the df form the full model and from the Ref-only model (let me add: I think it is a valid argument that the df from the latter should be used but I am sure we can have a good fight over that some other time).


Yep. 71.

❝ A=read.table("EMAset1b.csv", header=T)


Can you post your EMAset1b.csv? Since my crystal ball is at the laundry I can’t check your code.
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
d_labes
★★★

Berlin, Germany,
2013-04-09 11:01
(4418 d 14:43 ago)

@ ElMaestro
Posting: # 10379
Views: 5,928
 

 EMA dataset 1 - programmers trouble

Post reply
Hi ElMaestro,

❝ Here's some code bits:

❝

❝ A=read.table("EMAset1b.csv", header=T)

❝ ....

❝ lnU=lnPE+s*CritT*sqrt(0.5*(1/nBABA+1/nABAB)) ##can't get it right.

❝ lnL=lnPE-s*CritT*sqrt(0.5*(1/nBABA+1/nABAB)) ##can't get it right.


Try it with:
lnU=lnPE+s*CritT*sqrt(0.25*(1/nBABA+1/nABAB))
lnL=lnPE-s*CritT*sqrt(0.25*(1/nBABA+1/nABAB))

or even easier fit a model with intercept and get the confidence interval via call of confint(). No need to fiddle with n in sequences:
all <- lm(log(Data) ~ as.factor(formulation)+ as.factor(period) +as.factor(subject), data=A)
lnPE <- coef(all)[2]
CI   <- confint(all,level=0.9)[2,]

BTW: I have omitted sequence since it doesn't contribute to the mse to have such a term. Moreover I have all effects in lower letters (caps lock not arrested :-D).

My SAS gives (without back transformation):
                              Standard
  Parameter    Estimate          Error   t Value   Pr > |t|     90% Confidence Limits

  test-ref   0.14547372     0.04650869      3.13     0.0020     0.06864574   0.22230170

The latter R code gives:
> print(lnPE)
as.factor(tmt)T
      0.1454737
> print(CI)
       5 %       95 %
0.06864569 0.22230165
—
Regards,

Detlew
 
ElMaestro
★★★

Denmark,
2013-04-09 12:12
(4418 d 13:32 ago)

@ d_labes
Posting: # 10380
Views: 5,848
 

 Still not exactly sure

Post reply
Thank you, Helmut and d_labes, for your input.

I am not planning to via the in-built confidence interval function because I can port it directly to C when I know exactly what's going on.

Helmut, I downloaded your spreadsheet, saved the dataset as csv with tab delimiter and adjusted the code to fit your format.

Here's some code:
A=read.table("EMAset1c.csv", header=T)

Model1=lm(log(A$Data)~0+factor(A$Formulation)+
factor(A$Period)+factor(A$Sequence)+factor(A$Subject))

lnPE=(coef(Model1)[2]-coef(Model1)[1])


B=subset(A, Formulation=='R')
Model2=lm(log(B$Data)~factor(B$Period)+factor(B$Sequence)+factor(B$Subject))
s=summary(Model2)$sigma
CV=sqrt(exp(s*s)-1.0)
DF=Model2$df.residual
CritT=qt(1-0.05, DF)

nTRTR=0
nRTRT=0

for (i in min(B$Subject):max(B$Subject))
if (i %in% B$Subject) ##if the subject exists
{
  C=subset(B, Subject==i)
  r=subset(C, Formulation=='R')
  if(length(r$Formulation)==2) ##if the subject ate pocus two times.
   {
    if (C$Sequence[1]=='TRTR') nTRTR=nTRTR+1
    if (C$Sequence[1]=='RTRT') nRTRT=nRTRT+1
   }
} ## there's prolly a more efficient way of doing it.

lnU=lnPE+s*CritT*sqrt(0.25*(1/nTRTR+1/nRTRT))
lnL=lnPE-s*CritT*sqrt(0.25*(1/nTRTR+1/nRTRT))
cat("PE=", exp(lnPE),"\n")
cat("Upper Limit=", exp(lnU),"\n")
cat("Lower Limit=", exp(lnL),"\n")
cat("CV=", CV, "\n")
cat("nTRTR , nRTRT:", nTRTR, ",", nRTRT, "\n")
 

I am still not really there :-D
—
Pass or fail!
ElMaestro
 
d_labes
★★★

Berlin, Germany,
2013-04-09 17:57
(4418 d 07:47 ago)

@ ElMaestro
Posting: # 10383
Views: 5,816
 

 Still not exactly sure - why

Post reply
Dear ElMaestro,

could you kindly give a more elaborate description of what is your problem?
Without charging me with the burden to run your code and figure out what's going on!
—
Regards,

Detlew
 
ElMaestro
★★★

Denmark,
2013-04-09 18:39
(4418 d 07:05 ago)

@ d_labes
Posting: # 10384
Views: 5,718
 

 My problem(s)

Post reply
Hi d_labes,

❝ could you kindly give a more elaborate description of what is your problem?

❝ Without charging me with the burden to run your code and figure out what's going on!


I would like to reproduce the result of EMA's Q&A with method A: 107.11% - 124.89% for the CI using R.
The code was an attempt at that and it was of course a miserable failure which can probably be explained by my lack of skills. :vomit:
—
Pass or fail!
ElMaestro
 
ElMaestro
★★★

Denmark,
2013-04-10 01:29
(4418 d 00:15 ago)

@ d_labes
Posting: # 10386
Views: 5,795
 

 More confusion

Post reply
Gentlemen...

I have the point estimate and s right, this much is certain.
Since something is wrong, I just need to figure out which n-per-sequence that would result in limits such as the lower one of 107.11:

for (n1 in 30:80)
for (n2 in 30:80)
 {
  df=n1+n2-2
  ct = qt(1-0.05, df)
   lnL=lnPE-s*ct*sqrt(0.25*(1/n1+1/n2))
  if (exp(lnL)>1.0710)
  if (exp(lnL)<1.0712) cat("n1=", n1, "n2=",n2, "lower limit=", exp(lnL), "\n")
 }


The output is:
n1= 34 n2= 77 lower limit= 1.071095
n1= 35 n2= 72 lower limit= 1.071009
n1= 35 n2= 73 lower limit= 1.0712
n1= 36 n2= 69 lower limit= 1.071178
n1= 37 n2= 65 lower limit= 1.07102
n1= 39 n2= 60 lower limit= 1.071096
n1= 40 n2= 58 lower limit= 1.071152
n1= 41 n2= 56 lower limit= 1.071138
n1= 42 n2= 54 lower limit= 1.071052
n1= 44 n2= 51 lower limit= 1.071036
n1= 45 n2= 50 lower limit= 1.071145
n1= 50 n2= 45 lower limit= 1.071145
n1= 51 n2= 44 lower limit= 1.071036
n1= 54 n2= 42 lower limit= 1.071052
n1= 56 n2= 41 lower limit= 1.071138
n1= 58 n2= 40 lower limit= 1.071152
n1= 60 n2= 39 lower limit= 1.071096
n1= 65 n2= 37 lower limit= 1.07102
n1= 69 n2= 36 lower limit= 1.071178
n1= 72 n2= 35 lower limit= 1.071009
n1= 73 n2= 35 lower limit= 1.0712
n1= 77 n2= 34 lower limit= 1.071095

So, realistic n1 and n2 add up to irrelevant values, therefore I might be inclined to believe the problem must be one of these lines:
df=n1+n2-2
ct = qt(1-0.05, df)
lnL=lnPE-s*ct*sqrt(0.25*(1/n1+1/n2))

:confused::confused::confused:
—
Pass or fail!
ElMaestro
 
d_labes
★★★

Berlin, Germany,
2013-04-10 12:51
(4417 d 12:53 ago)

@ ElMaestro
Posting: # 10396
Views: 5,717
 

 Confusion help

Post reply
My Dear,

don't puzzle in one of your mentioned directions.
The formulas used are obviously only correct in case of no missings. And the EMA set I has missings.

In case of missings you have to use a complicated formula involving the design matrix only ElMaestro can understand :-D. Have a look at the explanation given for the ESTIMATE statement of SAS Proc GLM.
But don't ask me what the vectors and matrices are there. I myself could it only understand after some pills of Schützomycin.

Hope this helps against confusion seizures.
—
Regards,

Detlew
 
UA Flag
Activity
 Admin contact
23,424 posts in 4,927 threads, 1,670 registered users;
191 visitors (0 registered, 191 guests [including 12 identified bots]).
Forum time: 01:45 CEST (Europe/Vienna)

Only dead fish go with the current.    Scuba divers' proverb

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