Simulation of intra-subject variability [Regulatives / Guidelines]

posted by ElMaestro  – Denmark, 2011-11-08 12:02 (4926 d 13:10 ago) – Posting: # 7634
Views: 8,934

(edited on 2011-11-08 17:16)

All hands on deck!

I hope this illustrates the issue of no treament differences when fixed effects between subjects are zerod or some scatter in fixed effects between subjects exists:
Per=c(1,2,3, 1,2,3, 1,2,3, 1,2,3, 1,2,3, 1,2,3, 1,2,3, 1,2,3, 1,2,3, 1,2,3, 1,2,3, 1,2,3)
Subj=c(1,1,1, 2,2,2, 3,3,3, 4,4,4, 5,5,5, 6,6,6, 7,7,7, 8,8,8, 9,9,9, 10,10,10, 11,11,11, 12,12,12)
Trt=c(1,2,3,1,2,3,1,3,2,1,3,2,2,1,3,2,1,3,2,3,1,2,3,1,3,2,1,3,2,1,3,1,2,3,1,2)
Seq=c("ABC","ABC","ABC","ABC","ABC","ABC","ACB","ACB","ACB","ACB","ACB",
"ACB","BAC","BAC","BAC","BAC","BAC","BAC","BCA","BCA","BCA","BCA","BCA",
"BCA","CBA","CBA","CBA","CBA","CBA","CBA","CAB","CAB","CAB","CAB","CAB","CAB")
Seq=sort(rep(1:6, 6))

SubjE1=c(rep(0,12)) ##fixed effects for subjs: zero for all 12
SubjE2=c(9*runif(12)) ##variation in fixed effects between subjects

TrtE=c(9.876, 8.765, 7.654) ## some fixed treatment effects
PerE=c(1.5, 1.5, 1.5) ##some fixed period effects all same
SeqE=c(2.0,2.0, 2.0, 2.0, 2.0, 2.0) ##some fixed sequence effects all same

set.seed(12041958)
rany = rnorm(36, mean=0, sd=2) ##our single error
Y1=c(1:36)
Y2=c(1:36)
for (i in 1:36)
 {
  indexTrt=Trt[i];
  indexSubj=Subj[i];
  indexPer=Per[i];
  indexSeq=Seq[i];
  rany =  rnorm(1, mean=0, sd=2)
  Y1[i] = TrtE[indexTrt] + SeqE[indexSeq] + SubjE1[indexSubj] + PerE[indexPer]+rany
  Y2[i] = TrtE[indexTrt] + SeqE[indexSeq] + SubjE2[indexSubj] + PerE[indexPer]+rany
 }
m1=lm(Y1~0+as.factor(Trt) + as.factor(Seq)+ as.factor(Subj) + as.factor(Per))
m2=lm(Y2~0+as.factor(Trt) + as.factor(Seq)+ as.factor(Subj) + as.factor(Per))
anova(m1)
cat("m1 no subject effect. They do not differ in this scenario.\n")
anova(m2)
cat("m2 p significant for subjects. They differ in this scenario.\n")
cat("m1 models all subjects fixed effects as zero.\n")
cat("m1 difference between treat effects 1 and 3:", coef(m1)[1]-coef(m1)[3], "\n")
cat("m1 residual error:", summary(m1)$sigma, "\n")
cat("m3 models subjects effects as fixed but different between subjects.\n")
cat("m2 difference between treat effects 1 and 3:", coef(m2)[1]-coef(m2)[3], "\n")
cat("m2 residual error:", summary(m2)$sigma, "\n")

#
#


Have a good day. As always it is prolly wrong and bugged.

:pirate:

EM

Pass or fail!
ElMaestro

Complete thread:

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

My doctor gave me six months to live,
but when I couldn’t pay the bill
he gave me six months more.    Walter Matthau

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