1 kawaguchi p.1/81

Size: px

Start display at page:

Download "1 kawaguchi p.1/81"

きみおつつの
5 years ago
Views:

1 1 kawaguchi p.1/81

2 AUC 4.4 p.2/81

3 X Z X = α + βz + e α : Z = 0 X ( ) β : Z X ( ) e : 0 σ 2 p.3/81

4 2.1 Z X Z : X : Beckmann p.4/81

5 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

6 x z p.6/81

7 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

8 R ( ) lm(x z) Call: lm(formula = x z) Coefficients: (Intercept) z p.8/81

9 SAS ( ) data beckmann; input z datalines; run; proc reg data=beckmann; model x=z ; run; p.9/81

10 SAS ( ) SAS p.10/81

11 X = α + β(z Z) + e α β a = X, b = (Zi Z)(X i X) (Zi Z) 2 Z X a b p.11/81

12 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

13 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

14 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

15 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

16 (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

17 (2) H 0 : β = 0 ( α%) t = b se(b) ( ) t > t n 2,α/2 = H 0 reject (H 0 reject 0 ) p.17/81

18 t t = b se(b) = r n 2 1 r 2 H 0 : β = 0 H 0 : ρ = 0 ρ p.18/81

19 (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

20 (2) H 0 : α = 0 ( α%) t = a se(a ) t > t n 2,α/2 = H 0 reject p.20/81

21 R ( ) coefficients(summary(lm(x z))) Estimate Std.Error t value Pr(> t ) (Intercept) z e p.21/81

22 R ( ) confint(lm(x z)) 2.5 % 97.5 % (Intercept) z p.22/81

23 SAS ( ) (CLB) proc reg data=beckmann; model x=z / CLB; run; p.23/81

24 (1) X = βz + e β ˆβ = Zi X i Z 2 i p.24/81

25 (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

26 R ( ) summary(lm(x z-1)) confint(lm(x z-1)) Estimate Std.Error t value Pr(> t ) z < β p.26/81

27 SAS ( ) (NOINT) proc reg data=beckmann; model x=z /NOINT CLB; run; p.27/81

28 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

29 ( 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

30 ( 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

31 α + β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

32 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

33 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

34 2.8 x z p.34/81

35 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

36 x z p.36/81

37 ( ) ( 2.1 ) X( ) Z( ) Ẑ = Z + X X b Z = Fieller p.37/81

38 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

39 Fieller θ = µ/ν E[U θv ] = 0, ˆV [U θv ] = s 2 (a θ 2 a 2 2) θ 100(1 α)% { (U θv ) 2 } θ : s 2 (a θ2 a 2 2 ) t2 m,α/2 t m,α/2 s 2 a 2 2 /V 2 < 1 p.39/81

40 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

41 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

42 R (Z ) new.x< 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

43 Z ( ) Z > hat.z Z 95% > m1+m2*sqrt(m31+m32) > m1-m2*sqrt(m31+m32) Z p.43/81

44 AUC AUC(Area Under the bood - level Curve) concentrations concentrations concentrations time time time Subject 1 Subject 2 Subject 3 p.44/81

45 SAS NLMIXED Example (Pinheiro and Bates, 1995) time( ) conc( ) dose( ) AUC R PK AUC p.45/81

46 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

47 log(auc) = α + β log( ) + ε AUC = e α ( ) β e ε { β 1 = AUC β = 1 = AUC p.47/81

48 H 0 : β = 1 log(auc) log( ) log(auc) log( ) = α + β log( ) + ε β H 0 : β = 0 p.48/81

49 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

50 log(auc) log(dose) log(auc) = log( ) p.50/81

51 R (H 0 : β = 1 ) coefficients(summary( lm((log(auc$inf)-log(dose)) log(dose)))) Estimate Std. Error t value Pr(> t ) (Intercept) <0.001 log(dose) p.51/81

52 R (β ) confint(lm((log(auc$inf)) log(dose))) 2.5 % 97.5 % (Intercept) log(dose) p.52/81

53 R ( ) plot(lm(log(auc$inf) log(dose))) Q-Q Cook p.53/81

54 Scale-Location plot Standardized residuals Fitted values 1 Cook s distance plot Obs. number p.54/81 Cook s distance

55 ( ) log(auc) log(dose) log(auc) = log( ) p.55/81

56 H 0 : β = 1 ( ) coefficients(summary( lm((log(auc$inf)[-1]-log(dose)[-1]) log(dose)[-1]))) Estimate Std. Error t value Pr(> t ) (Intercept) <0.001 log(dose)[-1] p.56/81

57 β ( ) confint( lm((log(auc$inf)[-1]) log(dose)[-1])) 2.5 % 97.5 % (Intercept) log(dose)[-1] p.57/81

58 log(auc) = log( ) H 0 : β = β 1.52 AUC (Subject 1 ) log(auc) = log( ) H 0 : β = β 1.44 AUC p.58/81

59 (Subject1 ) Scale-Location plot Cook s distance plot Standardized residuals Cook s distance Fitted values Obs. number p.59/81

60 X: Z: D Z i = log 2 ( /3.5) (Z i ) n i Xi s i 1 3.5(0) (1) (2) s i : X p.60/81

61 Z 2 0 Z 0 p.61/81

62 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 p.62/81

63 X Z X = α + βz + γz 2 + ε ( ) s s 2 = (ni 1)s 2 i (ni 1) p.63/81

64 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

65 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

66 R ( : ) > a [,1] [1,] 1.80 [2,] 2.21 [3,] a = 1.80 b = 2.21 c = 0.70 p.66/81

67 2 H 0 : Z 2 γ = 0 (α/2%) L c = c se(c) se(c) = s m (3,3) m (3,3) M L c > t n 3,α/4 = H 0 reject p.67/81

68 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, )# > p.c H 0 : γ = % p.68/81

69 (1) X Z : Z a X Z = a + bz + cz 2 = Z b c Z = (1,Z,Z 2 ) p.69/81

70 (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

71 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

72 x z p.72/81

73 X = α + βz + ε ˆβ = ni (Z i Z)(X i X) ni (Z i Z) 2 ˆα = ni Xi ˆβ n i Z i ni p.73/81

74 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] > alpha [1] ˆβ = 0.786, ˆα = p.74/81

75 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

76 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, ) p.76/81

77 : > L.beta > p.beta H 0 : β = % p.77/81

78 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

79 x z p.79/81

80 4.6 Z 2 0 Z 0 Z D p.80/81

81 J.L. (KR ) (2004). 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, p.81/81

(lm) lm AIC 2 / 1

W707 s-taiji@is.titech.ac.jp 1 / 1 (lm) lm AIC 2 / 1 : y = β 1 x 1 + β 2 x 2 + + β d x d + β d+1 + ϵ (ϵ N(0, σ 2 )) y R: x R d : β i (i = 1,..., d):, β d+1 : ( ) (d = 1) y = β 1 x 1 + β 2 + ϵ (d > 1) y