1 kawaguchi atsushi@kurume-u.ac.jp 2005 7 2 p.1/81
2.1 2.2 2.2.3 2.3 AUC 4.4 p.2/81
X Z X = α + βz + e α : Z = 0 X ( ) β : Z X ( ) e : 0 σ 2 p.3/81
2.1 Z X 1 0.045 2 0.114 4 0.215 6 0.346 7 0.41 8 0.52 Z : X : Beckmann 10 0.67 15 0.942 p.4/81
R ( ) z<-c(1, 2, 4, 6, 7, 8, 10, 15) x<-c(0.045, 0.114, 0.215, 0.346, 0.410, 0.520, 0.670, 0.942) z, x plot(z,x) p.5/81
x 0.2 0.4 0.6 0.8 2 4 6 8 10 12 14 z p.6/81
n (Z 1,X 1 ),...,(Z n,x n ) S 2 = n (X i α βz i ) 2 i=1 a b b = (Zi Z)(X i X) (Zi Z) 2 a = X b Z p.7/81
R ( ) lm(x z) Call: lm(formula = x z) Coefficients: (Intercept) z -0.02777 0.06574 p.8/81
SAS ( ) data beckmann; input z x @@; datalines; 1 0.045 2 0.114 4 0.215 6 0.346 7 0.410 8 0.520 10 0.670 15 0.942 run; proc reg data=beckmann; model x=z ; run; p.9/81
SAS ( ) SAS p.10/81
X = α + β(z Z) + e α β a = X, b = (Zi Z)(X i X) (Zi Z) 2 Z X a b p.11/81
s ZX Z X s ZX = (Zi Z)(X i X) n 1 b = s ZX s X Z X r ( 1 r 1) r = s ZX s Z s X b = r s X s Z s Z = 1 n 1 (Zi Z) 2 s X = 1 n 1 (Xi X) 2 p.12/81
R 2 = S2 1 S 2 2 S 2 1 S 2 1 = (X i X) 2, S 2 2 = [X i X b(z i Z)] 2 R 2 = r 2 p.13/81
b β U,T R Z X = U +f, Z = T +e, e f R Z = σ2 T σ 2 T + σ2 e, R X = σ2 U σ 2 U + σ2 f β U,T U V Z Z p.14/81
e σ 2 s 2 = = (Xi a bz i ) 2 n 2 (Xi X a b(z i Z)) 2 n 2 (E[s 2 ] = σ 2 ) n 2 p.15/81
(1) β 100(1 α)% b t n 2,α/2 se(b) β b + t n 2,α/2 se(b) se(b) = s (Zi Z) 2 s t ν,p ν t p% p.16/81
(2) H 0 : β = 0 ( α%) t = b se(b) ( ) t > t n 2,α/2 = H 0 reject (H 0 reject 0 ) p.17/81
t t = b se(b) = r n 2 1 r 2 H 0 : β = 0 H 0 : ρ = 0 ρ p.18/81
(1) α 100(1 α)% a t n 2,α/2 se(a ) α a + t n 2,α/2 se(a ) se(a ) = s 1 n + Z2 (Zi Z) 2 p.19/81
(2) H 0 : α = 0 ( α%) t = a se(a ) t > t n 2,α/2 = H 0 reject p.20/81
R ( ) coefficients(summary(lm(x z))) Estimate Std.Error t value Pr(> t ) (Intercept) -0.028 0.017-1.668 0.146 z 0.066 0.002 31.063 7.39e-08 0 0 p.21/81
R ( ) confint(lm(x z)) 2.5 % 97.5 % (Intercept) -0.069 0.013 z 0.061 0.071 0 0 p.22/81
SAS ( ) (CLB) proc reg data=beckmann; model x=z / CLB; run; p.23/81
(1) X = βz + e β ˆβ = Zi X i Z 2 i p.24/81
(2) β 100(1 α)% ˆβ t n 2,α/2 se(ˆβ) β ˆβ + t n 2,α/2 se(ˆβ) se(β) = v Z 2 i v = 1 n 1 (Xi ˆβZ i ) 2 p.25/81
R ( ) summary(lm(x z-1)) confint(lm(x z-1)) Estimate Std.Error t value Pr(> t ) z 0.063 0.001 49.11 <0.001 0.060 β 0.066 p.26/81
SAS ( ) (NOINT) proc reg data=beckmann; model x=z /NOINT CLB; run; p.27/81
Z = Z 0 α + βz 0 100(1 α)% a + bz 0 ± t n 2,α/2 s 1 n + (Z 0 Z) 2 (Zi Z) 2 Z 0 = 0 p.28/81
( 1) Z = Z i (i = 1, 2,... ) = (Bonferroni ) P {a + bz (α + βz)} ( ) 2 M s 2 1 (Z Z) 2 α, Z = 1 α n + (Zi Z) 2 M α p.29/81
( 2) P max Z {a + bz (α + βz)} ( ) 2 s 2 1 n + (Z Z) 2 (Zi Z) 2 M α = 1 α M α = ( X E[ X]) 2 σ 2 /n + (b β)2 σ 2 / (Z i Z) 2 s 2 /σ 2 2F 2,n 2 p.30/81
α + βz 100(1 α)% L(Z ) < α + βz < U(Z ) Z Z L(Z ) = X + b(z Z) A(Z ) U(Z ) = X + b(z Z) + A(Z ) ( 1 A(Z ) = 2F 2,n 2,α s 2 n + (Z Z) ) 2 (Zi Z) 2 p.31/81
R ( 2.8) 1 n<-length(z) new.z<- seq(2, 12, by=0.02) ones<-rep(1, length(new.z)) pre.z<-rbind(ones, new.z) a<-lm(x z)$coefficients plot(new.z, a%*%pre.z, ylim=c(0,1), type="l", lwd=2, xlab="z", ylab="x") p.32/81
R ( 2.8) 2 sz<-sum((z-mean(z))ˆ2) s2<-sum((lm(x z)$residuals)ˆ2)/(n-2) s1<-1/n+(new.z-mean(z))ˆ2/sz U<-a%*%pre.z+sqrt(2*qf(0.95,2,6)*s2*s1) L<-a%*%pre.z-sqrt(2*qf(0.95,2,6)*s2*s1) points(new.z, U, type="l") points(new.z, L, type="l") points(z, x) p.33/81
2.8 x 0.0 0.2 0.4 0.6 0.8 1.0 2 4 6 8 10 12 z p.34/81
R ( ) p.u<-a%*%pre.z+sqrt(qt(0.95,2,6)*s2*s1) p.l<-a%*%pre.z-sqrt(qt(0.95,2,6)*s2*s1) points(new.z, p.u, type="l",lty=3) points(new.z, p.l, type="l",lty=3) p.35/81
x 0.0 0.2 0.4 0.6 0.8 1.0 2 4 6 8 10 12 z p.36/81
( ) ( 2.1 ) X( ) Z( ) Ẑ = Z + X X b Z = Fieller p.37/81
Fieller U N(µ,sa 1 ), V N(ν,sa 2 ), s : g = t2 m,α/2 s2 a 2 2 V 2 g < 1 (V 0 ) µ/ν 100(1 α)% ( 1 U 1 g V ± t ) m,α/2s (1 g)a 2 1 V + U2 a 2 2 V 2 p.38/81
Fieller θ = µ/ν E[U θv ] = 0, ˆV [U θv ] = s 2 (a 2 1 + θ 2 a 2 2) θ 100(1 α)% { (U θv ) 2 } θ : s 2 (a 2 1 + θ2 a 2 2 ) t2 m,α/2 t m,α/2 s 2 a 2 2 /V 2 < 1 p.39/81
Z (Fieller ) U = X X V = b n+1 s se(u) = s n se(v ) = (Zi Z) 2 = a 2 1 = (n + 1)/n a 2 2 = 1/ (Z i Z) 2 g = t2 m,α/2 s2 a 2 2 V 2 = t2 n 2,α/2 s2 b 2 (Z i Z) 2 p.40/81
Z (Fieller ) Fieller = (Z Z) 100(1 α)% E[U] E[V ] = β(z Z) β X X b(1 g) ±t n 2,α/2s b(1 g) (1 g)(n + 1) n + (X X) 2 b 2 (Z i Z) 2 Z Z 100(1 α)% p.41/81
R (Z ) new.x<-0.577 hat.z<-mean(z)+(new.x-mean(x))/a[2] g1<-s2*(qt(0.975, (length(z)-2)))ˆ2 g2<-a[2]ˆ2*sum((z-mean(z))ˆ2) g<-g1/g2 m1<-mean(z)+(new.x-mean(x))/(a[2]*(1-g)) m2<-qt(0.975, (n-2))*sqrt(s2)/(a[2]*(1-g)) m31<-((1-g)*(n+1))/(n) m32<-(new.x-mean(x))ˆ2/(a[2]ˆ2*sz) m1+m2*sqrt(m31+m32) m1-m2*sqrt(m31+m32) p.42/81
Z ( ) Z > hat.z 9.199561 Z 95% > m1+m2*sqrt(m31+m32) 10.24142 > m1-m2*sqrt(m31+m32) 8.189849 8.19 Z 10.24 p.43/81
AUC AUC(Area Under the bood - level Curve) concentrations 0 2 4 6 8 10 12 concentrations 0 2 4 6 8 10 12 concentrations 0 2 4 6 8 10 12 0 5 10 15 20 25 time 0 5 10 15 20 25 time 0 5 10 15 20 25 time Subject 1 Subject 2 Subject 3 p.44/81
SAS NLMIXED Example (Pinheiro and Bates, 1995) time( ) conc( ) dose( ) 12 25 11 AUC R PK AUC p.45/81
R (AUC ) data<-read.table("auc_data.txt", header=true) library(pk) auc<-matrix(0, 12, 3) for(i in 1:12) {auc[i,]<- AUC(conc=data$CONC[data$SUBJECT==i], time=data$time[data$subject==i])$auc} colnames(auc) <- c("obs","inter","inf") auc<-data.frame(auc) p.46/81
log(auc) = α + β log( ) + ε AUC = e α ( ) β e ε { β 1 = AUC β = 1 = AUC p.47/81
H 0 : β = 1 log(auc) log( ) log(auc) log( ) = α + β log( ) + ε β H 0 : β = 0 p.48/81
R ( ) dose<-numeric(0) for(i in 1:12) {dose[i]<-data$dose[data$subject==i]} ones<-rep(1, length(dose)) predictor<-rbind(ones, log(dose)) plot(log(dose), log(auc$inf), xlab="log(dose)",ylab="log(auc)",ylim=c(4.4, 5.4)) points(log(dose), lm(log(auc$inf) log(dose))$coe%*%predictor,type="l") p.49/81
log(auc) 4.4 4.6 4.8 5.0 5.2 5.4 1.2 1.3 1.4 1.5 1.6 1.7 log(dose) log(auc) = 4.13 + 0.42 log( ) p.50/81
R (H 0 : β = 1 ) coefficients(summary( lm((log(auc$inf)-log(dose)) log(dose)))) Estimate Std. Error t value Pr(> t ) (Intercept) 4.13 0.75 5.47 <0.001 log(dose) 0.58 0.49 1.17 0.27 1 p.51/81
R (β ) confint(lm((log(auc$inf)) log(dose))) 2.5 % 97.5 % (Intercept) 2.44 5.81 log(dose) 0.68 1.52 p.52/81
R ( ) plot(lm(log(auc$inf) log(dose))) Q-Q Cook p.53/81
Scale-Location plot 1 11 10 Standardized residuals 0.0 0.5 1.0 1.5 4.60 4.70 4.80 Fitted values 1 Cook s distance plot 10 6 2 4 6 8 10 12 Obs. number p.54/81 Cook s distance 0.0 0.1 0.2 0.3 0.4 0.5
( ) log(auc) 4.4 4.6 4.8 5.0 5.2 5.4 1.2 1.3 1.4 1.5 1.6 1.7 log(dose) log(auc) = 3.59 + 0.73 log( ) p.55/81
H 0 : β = 1 ( ) coefficients(summary( lm((log(auc$inf)[-1]-log(dose)[-1]) log(dose)[-1]))) Estimate Std. Error t value Pr(> t ) (Intercept) 3.59 0.48 7.44 <0.001 log(dose)[-1] 0.27 0.31 0.86 0.41 1 p.56/81
β ( ) confint( lm((log(auc$inf)[-1]) log(dose)[-1])) 2.5 % 97.5 % (Intercept) 2.50 4.69 log(dose)[-1] 0.02 1.44 p.57/81
log(auc) = 4.13 + 0.42 log( ) H 0 : β = 1 0.68 β 1.52 AUC (Subject 1 ) log(auc) = 3.59 + 0.73 log( ) H 0 : β = 1 0.02 β 1.44 AUC p.58/81
(Subject1 ) Scale-Location plot Cook s distance plot 9 8 10 8 9 10 Standardized residuals 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Cook s distance 0.0 0.5 1.0 1.5 2.0 2.5 4.5 4.6 4.7 4.8 4.9 2 4 6 8 10 Fitted values Obs. number p.59/81
4.4 4.6 X: Z: D Z i = log 2 ( /3.5) (Z i ) n i Xi s i 1 3.5(0) 10 1.8 1.09 2 7.0(1) 8 3.31 1.63 3 14.0(2) 13 3.42 1.02 s i : X p.60/81
Z 2 0 Z 0 p.61/81
X i : Z i X (i = 1, 2,...,g) H 0 : 2 F = ni ( X i X i) 2 (g 3)s 2 F > F g 3,ν,α/3 = H 0 reject 4.6 3 p.62/81
X Z X = α + βz + γz 2 + ε ( ) s s 2 = (ni 1)s 2 i (ni 1) p.63/81
a b c a ni Xi b = M 1 ni Xi Z i c ni Xi Zi 2 ni ni Z i ni Zi 2 M = ni Z i ni Zi 2 ni Zi 3 ni Zi 2 ni Zi 3 ni Zi 4 p.64/81
R ( ) z<-c(0, 1, 2) n<-c(10, 8, 13) x<-c(1.80, 3.31, 3.42) M<-c(sum(n), sum(n*z), sum(n*zˆ2), sum(n*z), sum(n*zˆ2), sum(n*zˆ3), sum(n*zˆ2), sum(n*zˆ3), sum(n*zˆ4)) M<-matrix(M, ncol=3) M2<-c(sum(n*x), sum(n*x*z), sum(n*x*zˆ2)) a<-solve(m)%*%m2 p.65/81
R ( : ) > a [,1] [1,] 1.80 [2,] 2.21 [3,] -0.70 a = 1.80 b = 2.21 c = 0.70 p.66/81
2 H 0 : Z 2 γ = 0 (α/2%) L c = c se(c) se(c) = s m (3,3) m (3,3) M 1 3 3 L c > t n 3,α/4 = H 0 reject p.67/81
R (2 ) s2<-((n-1)%*%sˆ2)/sum(n-1) se.c<-sqrt(s2*solve(m)[3,3]) L.c<-a[3]/se.c p.c<-1-pt(abs(l.c), 28, 0.0125)# > p.c 0.089 H 0 : γ = 0 0.025% p.68/81
(1) X Z : Z a X Z = a + bz + cz 2 = Z b c Z = (1,Z,Z 2 ) p.69/81
(2) 2 Working-Hotelling ( 100(1 α)%) X Z ± 3F 3,n. g,α se(x z) n. = n i se(x Z) = s Z M 1 Z p.70/81
R ( ) new.z<-seq(0, 3, 0.2) ones<-rep(1, length(new.z)) new.z<-cbind(ones, new.z, new.zˆ2) se.x<-sqrt(s2*diag(new.z%*%solve(m)%*%t(new.z))) plot(new.z, new.z%*%a, ylim=c(0,6), type="l", xlab="z", ylab="x") U<-new.Z%*%a+sqrt(3*qf(0.95, 3, sum(n)-3))*se.x L<-new.Z%*%a-sqrt(3*qf(0.95, 3, sum(n)-3))*se.x points(new.z, U, type="l") points(new.z, L, type="l") points(z, x) p.71/81
x 0 1 2 3 4 5 6 0.0 0.5 1.0 1.5 2.0 2.5 3.0 z p.72/81
X = α + βz + ε ˆβ = ni (Z i Z)(X i X) ni (Z i Z) 2 ˆα = ni Xi ˆβ n i Z i ni p.73/81
beta<-((sum(n)*sum(n*x*z)-(sum(n*x))*(sum(n*z))) /(sum(n)*sum(n*zˆ2)-sum(n*z)ˆ2)) alpha<-((sum(n*x)-beta*sum(n*z))/(sum(n))) > beta [1] 0.786 > alpha [1] 2.007 ˆβ = 0.786, ˆα = 2.007 p.74/81
H 0 : β = 0 ( α/2%) Lˆβ = ˆβ se(ˆβ) se(ˆβ) = s ni ni (Z i Z) 2 Lˆβ > t n 3,α/4 = H 0 reject p.75/81
R ( ) se.beta<-(sqrt((s2*sum(n))/ (sum(n)*sum(n*zˆ2)-sum(n*z)ˆ2))) L.beta<-beta/se.beta p.beta<-1-pt(l.beta, 28, 0.0125) p.76/81
: > L.beta 3.067 > p.beta 0.002 H 0 : β = 0 0.025% p.77/81
R ( ) a1<-coefficients(lm(x z, weight=n)) ones<-rep(1, length(new.z)) new.z1<-cbind(ones, new.z) points(new.z, new.z1%*%a1, type="l") p.78/81
x 0 1 2 3 4 5 6 0.0 0.5 1.0 1.5 2.0 2.5 3.0 z p.79/81
4.6 Z 2 0 Z 0 Z D p.80/81
J.L. (KR ) (2004). RjpWiki http://www.okada.jp.org/rwiki/index.php?rjpwiki Pinheiro, J.C. and Bates, D.M. (1995), "Approximations to the Log-likelihood Function in the Nonlinear Mixed-effects Model," Journal of Computational and Graphical Statistics, 4, 12-35. p.81/81