薬物動態の統計数理　医学統計ライブラリー　後藤昌司・大門貴志先生

```{r}

#コンパートメント・モデル　p.35

#4次のルンゲクッタ　連立微分方程式

k12=0.2 #薬物投与部位から体循環へのコンパートメントの移行速度定数
k20=0.5　#体循環コンパートメントからの移行速度定数

#dx1/dt=-k12*x1

f=function(tt,xx,yy){

return(-k12*xx)  

}

#dx2/dt=k12*x1-k20*x2

g=function(tt,xx,yy){

return(k12*xx-k20*yy)  

}

t=0
x=4.5 #初期の薬物投与量：これから投与する量
y=0　#初期の薬物循環量：0　

x_vec=c()
y_vec=c()
t_vec=c()

ite=2000
h=0.01

for(l in 1:ite){

x_vec=c(x_vec,x)
y_vec=c(y_vec,y)
t_vec=c(t_vec,t)  

s11=f(t,x,y)
s12=g(t,x,y)
s21=f(t+h/2,x+h*s11/2,y+h*s11/2)
s22=g(t+h/2,x+h*s12/2,y+h*s12/2)
s31=f(t+h/2,x+h*s21/2,y+h*s21/2)
s32=g(t+h/2,x+h*s22/2,y+h*s22/2)
s41=f(t+h,x+h*s31,y+h*s31)
s42=g(t+h,x+h*s32,y+h*s32)
x=x+h*(s11+2*s21+2*s31+s41)/6
y=y+h*(s12+2*s22+2*s32+s42)/6
t=t+h

}

#投与部位コンパートメントにおける時刻tでの薬物量
plot(t_vec,x_vec,xlab="t",ylab="x",xlim=c(min(t_vec),max(t_vec)),ylim=c(min(x_vec),max(x_vec)),col=2,main="投与部位コンパートメントにおける時刻tでの薬物量")

#体循環コンパートメントにおける時点tでの薬物量
#吸収=>代謝=>排泄のようす
plot(t_vec,y_vec,xlab="t",ylab="y",xlim=c(min(t_vec),max(t_vec)),ylim=c(min(y_vec),max(y_vec)),col=3,main="体循環コンパートメントにおける時点tでの薬物量")



```


```{r}

#薬物動態の統計数理　p.50

data(Theoph)

X=cbind(rep(1,nrow(Theoph)),Theoph$Dose)

Y=Theoph$conc
Y[Y==0]=0.0001

t=Theoph$Time

n=length(t)

p=1
b0=rep(0.01,p)
b1=rep(0.01,p)
lam=1
sig=sum(((Y^lam-1)/lam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)/n


ite=100000
eta=10^(-4)
eta_lam=10^(-4)
b0_pre=b0+1
b1_pre=b1+1
lam_pre=lam+1
Ylam=ifelse(rep(lam,length(Y))!=0,(Y^lam-1)/lam,log(Y))
loglik=-n*log(2*pi*sig)/2+(lam-1)*sum(log(Y))-sum((Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)/(2*sig)
loglik_pre=loglik+1

for(l in 1:ite){
  
Ylam=ifelse(rep(lam,length(Y))!=0,(Y^lam-1)/lam,log(Y))
loglik_pre=loglik+1
while(abs(loglik_pre-loglik)>0.00001)({  
loglik_pre=loglik
db0=t(exp(t(t(c(t)))%*%t(c(b1))))%*%(Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)/sig  
db1=b0*(t(exp(t(t(c(t)))%*%t(c(b1))))%*%(t*(Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)))/sig
b0=b0+eta*db0
b1=b1+eta*db1
sig=sum((Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)/n
loglik=-n*log(2*pi*sig)/2+(lam-1)*sum(log(Y))-sum((Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)/(2*sig)
print(loglik)
})  

print(loglik)
loglik_pre=loglik+1
while(abs(loglik_pre-loglik)>0.00001)({
loglik_pre=loglik
Ylam=ifelse(rep(lam,length(Y))!=0,(Y^lam-1)/lam,log(Y))
dlam=sum(log(Y))-sum((1/lam)*((Y^lam)*log(Y)-(Y^lam-1)/lam)*(Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0))/(sig)
lam=lam+eta_lam*dlam
sig=sum((Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)/n
loglik=-n*log(2*pi*sig)/2+(lam-1)*sum(log(Y))-sum((Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)/(2*sig)
print(loglik)
})

print(loglik)

}









```

```{r}

#薬物動態の統計数理　p.50

data(Theoph)

X=cbind(rep(1,nrow(Theoph)),Theoph$Dose)

Y=Theoph$conc
Y[Y==0]=0.0001

t=cbind(Theoph$Time,Theoph$Time^2)

n=length(t)

p=2
b0=rep(0.01,p)
b1=rep(0.01,p)
lam=1
if(p==1){
sig=sum(((Y^lam-1)/lam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)/n
}else{
sig=sum(((Y^lam-1)/lam-exp(t(t(t)*c(b1)))%*%b0)^2)/n
}

ite=100000
eta=10^(-6)
eta_lam=10^(-6)
b0_pre=b0+1
b1_pre=b1+1
lam_pre=lam+1
Ylam=ifelse(rep(lam,length(Y))!=0,(Y^lam-1)/lam,log(Y))
if(p==1){
Se=sum((Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)    
}else{
Se=sum((Ylam-exp(t(t(t)*c(b1)))%*%b0)^2)    
}

loglik=-n*log(2*pi*sig)/2+(lam-1)*sum(log(Y))-Se/(2*sig)
loglik_pre=loglik+1

for(l in 1:ite){
  
Ylam=ifelse(rep(lam,length(Y))!=0,(Y^lam-1)/lam,log(Y))
loglik_pre=loglik+1
while(abs(loglik_pre-loglik)>0.00001)({  
loglik_pre=loglik
if(p==1){
db0=t(exp(t(t(c(t)))%*%t(c(b1))))%*%(Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)/sig  
db1=b0*(t(exp(t(t(c(t)))%*%t(c(b1))))%*%(t*(Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)))/sig
}else{
db0=(exp(t(t)*c(b1)))%*%(Ylam-(exp(t(t(t)*c(b1))))%*%b0)/sig  
db1=b0*apply(t(exp(t(t)*c(b1)))*(t*c(Ylam-exp(t(t(t)*c(b1)))%*%b0)),2,sum)/sig
}
b0=b0+eta*db0
b1=b1+eta*db1
if(p==1){
sig=sum(((Y^lam-1)/lam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)/n
Se=sum((Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)
}else{
sig=sum(((Y^lam-1)/lam-(exp(t(t(t)*c(b1))))%*%b0)^2)/n
Se=sum((Ylam-(exp(t(t(t)*c(b1))))%*%b0)^2)   
}
loglik=-n*log(2*pi*sig)/2+(lam-1)*sum(log(Y))-Se/(2*sig)
#print(loglik)
})  

print(loglik)
loglik_pre=loglik+1
while(abs(loglik_pre-loglik)>0.00001)({
loglik_pre=loglik
Ylam=ifelse(rep(lam,length(Y))!=0,(Y^lam-1)/lam,log(Y))
if(p==1){
dlam=sum(log(Y))-sum((1/lam)*((Y^lam)*log(Y)-(Y^lam-1)/lam)*(Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0))/(sig)
}else{
dlam=sum(log(Y))-sum((1/lam)*((Y^lam)*log(Y)-(Y^lam-1)/lam)*(Ylam-(exp(t(t(t)*c(b1))))%*%b0))/(sig)
}
lam=lam+eta_lam*dlam
if(p==1){
sig=sum((Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)/n
Se=sum((Ylam-exp(t(t(c(t)))%*%t(c(b1)))%*%b0)^2)
}else{
sig=sum((Ylam-(exp(t(t(t)*c(b1))))%*%b0)^2)/n
Se=sum((Ylam-(exp(t(t(t)*c(b1))))%*%b0)^2)   
}
loglik=-n*log(2*pi*sig)/2+(lam-1)*sum(log(Y))-Se/(2*sig)
#print(loglik)
})

print(loglik)

}

```
```{r}

#薬物動態の統計数理　p.43

#べき乗型

data(Theoph)

X=cbind(rep(1,nrow(Theoph)),Theoph$Dose)

Y=Theoph$conc
Y[Y==0]=0.0001
t=cbind(Theoph$Time)
n=length(Y)



p=1
b0=rep(1,p)
b1=rep(0.001,p)
sita=rep(0.1,3)

ite=360000
eta=10^(-7)
for(l in 1:ite){
  
if(p==1){
fval=exp(t(t(c(t)))%*%t(c(b1)))%*%b0  
}else{
fval=exp(t(t(t)*c(b1)))%*%b0    
} 
sig=sum((Y-fval)^2/(n*(sita[1]+sita[2]*fval^sita[3])^2))
  
dfb0=t(exp(t(t)*b1))
dfb1=t(t(t)*exp(t(t)*b1)*b0)
dgsita=cbind(1,fval^(sita[3]),sita[2]*fval*log(fval))
dgb0=dfb0*c((sita[1]+sita[2]*fval^sita[3])*sita[3]/fval)  
dgb1=dfb1*c((sita[1]+sita[2]*fval^sita[3])*sita[3]/fval)

dlb0=-dgb0/c(sita[1]+sita[2]*fval^sita[3])+dgb0*c(((Y-fval)^2)/(sig*(sita[1]+sita[2]*fval^sita[3])^3))+dfb0*c((Y-fval)/((sita[1]+sita[2]*fval^sita[3])^2))
dlb0=apply(dlb0,2,sum)
dlb1=dgb1/c(sita[1]+sita[2]*fval^sita[3])+dgb1*c(((Y-fval)^2)/(sig*(sita[1]+sita[2]*fval^sita[3])^3))+dfb1*c((Y-fval)/((sita[1]+sita[2]*fval^sita[3])^2))
dlb1=apply(dlb1,2,sum)
dlsita=apply(-dgsita/c(sita[1]+sita[2]*fval^sita[3])+dgsita*c(((Y-fval)^2)/(sig*(sita[1]+sita[2]*fval^sita[3])^3)),2,sum)

b0=b0+eta*dlb0
b1=b1+eta*dlb1
sita=sita+eta*dlsita

loglik=-n*log(2*pi*sig)/2-sum(log((sita[1]+sita[2]*fval^sita[3])^2))/2-sum(((Y-fval)^2)/((sita[1]+sita[2]*fval^sita[3])^2))/(2*sig)

print(loglik)

}



```

```{r}

#薬物動態の統計数理　p.43

#指数型

data(Theoph)

X=cbind(rep(1,nrow(Theoph)),Theoph$Dose)

Y=Theoph$conc
Y[Y==0]=0.0001
t=cbind(Theoph$Time,Theoph$Time^2)
n=length(Y)



p=2
b0=rep(1,p)
b1=rep(0.001,p)
sita=rep(0.1,2)

ite=150000
eta=10^(-7)
for(l in 1:ite){
  
if(p==1){
fval=exp(t(t(c(t)))%*%t(c(b1)))%*%b0  
}else{
fval=exp(t(t(t)*c(b1)))%*%b0    
} 
sig=sum((Y-fval)^2/(n*(exp(sita[1]*fval^sita[2]))^2))
  
dfb0=t(exp(t(t)*b1))
dfb1=t(t(t)*exp(t(t)*b1)*b0)
dgsita=cbind(fval*exp(sita[1]*fval^sita[2]),(sita[1]*fval^sita[2])*log(fval)*exp(sita[1]*fval^sita[2]))
dgb0=prod(sita)*dfb0*c(exp(sita[1]*fval^sita[2])*fval^(sita[2]-1)) 
dgb1=prod(sita)*dfb1*c(exp(sita[1]*fval^sita[2])*fval^(sita[2]-1))

dlb0=-dgb0/c(exp(sita[1]*fval^sita[2]))+dgb0*c(((Y-fval)^2)/(sig*(exp(sita[1]*fval^sita[2]))^3))+dfb0*c((Y-fval)/((exp(sita[1]*fval^sita[2]))^2))
dlb0=apply(dlb0,2,sum)
dlb1=dgb1/c(exp(sita[1]*fval^sita[2]))+dgb1*c(((Y-fval)^2)/(sig*(exp(sita[1]*fval^sita[2]))^3))+dfb1*c((Y-fval)/((exp(sita[1]*fval^sita[2]))^2))
dlb1=apply(dlb1,2,sum)
dlsita=apply(-dgsita/c(exp(sita[1]*fval^sita[2]))+dgsita*c(((Y-fval)^2)/(sig*(exp(sita[1]*fval^sita[2]))^3)),2,sum)

b0=b0+eta*dlb0
b1=b1+eta*dlb1
sita=sita+eta*dlsita

loglik=-n*log(2*pi*sig)/2-sum(log((exp(sita[1]*fval^sita[2]))^2))/2-sum(((Y-fval)^2)/(exp(sita[1]*fval^sita[2])^2))/(2*sig)

print(loglik)

}



```