############################################################# # Complete link functions, output also includes first and # # second derivatives for Heissian matrix # ############################################################# ##### mm link ##### mmlink<-function(theta,j=1) { h=theta/(j+theta) h.1d=j/(theta+j)^2 h.2d=-2*j/(theta+j)^3 list(h=h,h.1d=h.1d,h.2d=h.2d) } ##### one parameter gamma link ##### ga1<-function(theta,j=1) { h=1/((1+j)^theta) h.1d=-1*log(1+j)/((1+j)^theta) h.2d=(log(1+j))^2/((1+j)^theta) list(h=h,h.1d=h.1d,h.2d=h.2d) } ##### Two-parameter gamma link ##### ga<-function(theta1,theta2=1, j){ h=1/(1+j*theta2)^theta1 h.1d=c(-1*h*log(1+j*theta2),-1*theta1*j*(1+theta2*j)^(-1*theta1-1)) h.2d.1=h*(log(1+j*theta2))^2 h.2d.4=theta1*(theta1+1)*(1+theta2*j)^(-1*theta1-2)*j^2 h.2d.2=j*(1+theta2*j)^(-1*theta1-1)*(theta1*log(1+theta2*j)-1) h.2d=matrix(c(h.2d.1,h.2d.2,h.2d.2,h.2d.4),2,2) list(h=h,h.1d=h.1d,h.2d=h.2d) } ##### Binomial link ##### bin<-function(theta,j=1) { h=theta^j h.1d=j*theta^(j-1) h.2d=j*(j-1)*theta^(j-2) list(h=h,h.1d=h.1d,h.2d=h.2d) } ##### Binomail power link ##### bin2<-function(theta,v,j=1){ h=theta^(j^v) h.1d=c(j^v*theta^(j^v-1), ifelse(j==0,0,h*log(theta)*log(j)*j^v)) h.2d.1=j^v*(j^v-1)*theta^(j^v-2) h.2d.4=ifelse(j==0,0,h.1d[2]*log(j)*(log(theta)*j^v+1)) h.2d.2=ifelse(j==0,0,h.1d[1]*log(j)*(1+log(theta))) h.2d=matrix(c(h.2d.1,h.2d.2,h.2d.2,h.2d.4),2,2) list(h=h, h.1d=h.1d, h.2d=h.2d) } ##### piecewise link ####### piecelink<-function(theta,j=1) { if (theta>1) { h=1/(1+j/2)^theta h.1d=-1*h*log(1+j/2) h.2d=h*(log(1+j/2))^2 } else { h= 2/(1+(1+j)^theta) h.1d=-1*h^2/2*(1+j)^theta*log(1+j) h.2d=h^2/2*(1+j)^theta*(log(1+j))^2*((1+j)^theta-1)/(1+(1+j)^theta) } list(h=h, h.1d=h.1d,h.2d=h.2d) } ##### powerstable link ####### pslink<-function(theta,j=1) { h=exp(-1*j^(1/(1+theta))) s2=j^(1/(1+theta)) h.1d=exp(-s2)*s2*log(j)/(1+theta)^2 h.2d=exp(-s2)*s2*log(j)*(-2*(1+theta)+(s2-1)*log(j))/(1+theta)^4 list(h=h,h.1d=h.1d,h.2d=h.2d) } #######link 6 ############# link6<-function(theta,j=1) { h=exp(-j^(1/(1+theta))) s2=j^(1/(1+theta)) h.1d=s2*log(j)/((1+s2)^2*(1+theta)^2) deri1=s2*log(j)/((1+s2)^2*(1+theta)^2) h.2d=-1*deri1*(2*(1+s2)*(1+theta)-(s2-1)*log(j))/((1+s2)*(1+theta)^2) list(h=h,h.1d=h.1d,h.2d=h.2d) } ########link 8 ###### link8<-function(theta,j=1) { h=(1+log(1+j))^(-theta) #delta=1 s2=1+log(1+j) h.1d=-s2^(-theta)*log(s2) h.2d=s2^(-theta)*(log(s2))^2 list(h=h,h.1d=h.1d,h.2d=h.2d) } #####link9 ############### link9<-function(theta,j=1) { h=(1+theta*log(1+j))^(-1) ##theta=1 s2=log(1+j) h.1d=-s2/(1+theta*s2)^2 h.2d=2*(s2)^2/(1+theta*s2)^3 list(h=h,h.1d=h.1d,h.2d=h.2d) } ###### link12 ############## link12<-function(theta,j=1) { h=1-(1-exp(-j))^(1/(1+theta)) es=exp(-j) etaInv=1/(1+theta) h.1d=(1-es)^etaInv*log(1-es)*etaInv^2 d1=(1-es)^etaInv*log(1-es)*etaInv^2 h.2d=-1*d1*(2+2*theta+log(1-es))*etaInv^2 list(h=h,h.1d=h.1d,h.2d=h.2d) } #####link from Qi and Xu (1998)paper########## qxlink<-function(theta,j=1) { h=ifelse(j==0, 1, 1/(1+theta)^j*(1-exp(-j))/j) h.1d=-1*(1-exp(-j))*(1+theta)^(-j-1) h.2d=(j+1)*(1-exp(-j))*(1+theta)^(-j-2) list(h=h,h.1d=h.1d,h.2d=h.2d) } ################################################################ # incomplete link functions, output includes first derivatives # # for calculating standard deviations of estimators # ################################################################ ##### incomplete gamma link ###### iga<-function(theta1, theta2, theta3, theta4,j=1) { h1=pgamma((theta4*j+1)*theta2,theta3)*gamma(theta3)-pgamma((theta4*j+1)*theta1,theta3)*gamma(theta3) h2=pgamma(theta2,theta3)*gamma(theta3)-pgamma(theta1,theta3)*gamma(theta3) tmp=(theta4*j+1)^theta3 if (h2==0){ h=theta1^j h.1d=c(binlink(theta1,j=j)$h.1d,0,0,0) } else { h=h1/(h2*tmp) h1d.1=1/tmp*1/h2^2*(dgamma((theta4*j+1)*theta1,theta3)*(-1)*(theta4*j+1)*h2-h1*(-1)*dgamma(theta1,theta3)) h1d.2=1/tmp*1/h2^2*(dgamma((theta4*j+1)*theta2,theta3)*(theta4*j+1)*h2-h1*dgamma(theta2,theta3)) h1d.3=h1/h2*ga(theta3,theta4,s=j)$h.1d[1]+1/tmp*1/h2^2*((gammapr(A=(theta4*j+1)*theta2, theta=theta3)$value-gammapr(A=(theta4*j+1)*theta1, theta=theta3)$value)*h2-h1*gammapr(A=theta2, theta=theta3)$value-gammapr(A=theta1, theta=theta3)$value) h1d.4=h1/h2*ga(theta3,theta4,s=j)$h.1d[2]+1/(tmp*h2)*(dgamma((theta4*j+1)*theta2,theta3)*j*theta2-dgamma((theta4*j+1)*theta1,theta3)*j*theta1) h.1d=c(h1d.1,h1d.2,h1d.3,h1d.4) } list(h=h, h.1d=h.1d) } ##### first derivative for gamma(theta) needed in iga ##### gammapr<-function(A=0,a2=Inf,theta) { gp<-integrate(function(x,theta1=theta) exp(-x)*x^(theta1-1)*log(x),A,a2)[1] return(gp) } ##### incomplete beta link ###### ibt<-function(theta1=1,theta2=0, theta3, theta4, j=1) { tmp1=beta(theta3+j, theta4) tmp2=beta(theta3, theta4) h1=pbeta(theta2, theta3+j, theta4)-pbeta(theta1, theta3+j, theta4) h2=pbeta(theta2, theta3, theta4)-pbeta(theta1, theta3, theta4) if (h2==0) { h=theta1^j h.1d=c(bin(theta1,s=j)$h.1d,0,0,0) } else { h=tmp1/tmp2*h1/h2 h.1d.1=tmp1/tmp2*1/h2^2*((-1)*dbeta(theta1, theta3+j, theta4)*h2-h1*(-1)*dbeta(theta1, theta3, theta4)) h.1d.2=tmp1/tmp2*1/h2^2*(dbeta(theta2, theta3+j, theta4)*h2-h1*dbeta(theta2, theta3, theta4)) h.1d.34=h1/(h2*tmp2^2)*(betapr(a1=(theta3+j),a2=theta4)*tmp2-tmp1*betapr(a1=theta3,a2=theta4))+tmp1/(tmp2*h2^2)*(ibtpr(a1= theta1,a2=theta2, b1=theta3+j, b2=theta4)*h2-h1*ibtpr(a1=theta1, a2=theta2, b1=theta3, b2=theta4)) h.1d=c(h.1d.1,h.1d.2,h.1d.34) } list(h=h, h.1d=h.1d) } ##### compute derivative of beta(a1, a2), the output is a vector ##### betapr<-function(a1,a2){ g1=gamma(a1) g2=gamma(a2) g3=gamma(a1+a2) d1=gammapr(theta=a1)$value d2=gammapr(theta=a2)$value d3=gammapr(theta=(a1+a2))$value pr1=1/(a1*a2)^2*(g1*g2*d3-g2*d1) pr2=1/(a1*a2)^2*(g1*g2*d3-g1*d2) dbt=c(pr1,pr2) return(dbt) } ##### derivatives for b1 and b2 in ibt w.r.t. integral limits ##### ibtpr<-function(a1,a2,b1,b2){ pr1<-integrate(function(p,t=b1) log(p)*p^(t-1)*(1-p)^(b2-1), a1,a2)[1] pr2<-integrate(function(p,t=b2) log(1-p)*p^(b1-1)*(1-p)^(t-1), a1,a2)[1] dbt=c(pr1$value,pr2$value) return(dbt) } ##### incomplete A- link ##### inA<-function(theta1, theta2, t=1) { tmp=log(theta2)-log(theta1) if (t==0) { h=1 h.1d=c(0,0) } else { h=(theta2^t-theta1^t)/(t*tmp) h.1d=c(((-1)*t*theta1^t*tmp+theta2^t-theta1^t)/(theta1*t*tmp^2),(t*theta2^t*tmp-theta2^t+theta1^t)/(theta2*t*tmp^2)) } list(h=h, h.1d=h.1d) } ######### loglikelihood function ########## loglike<-function(beta,data=E2,ep=1e-15,link=mmlink){ #returns (-1) log likehood of ET #data=(dose lts malf) N<- dim(data)[1] m<- data[,2] y<- data[,1] loglik<- 0 for(i in 1:N){ myi= m[i]-y[i] si<-0 for (k in 0:myi){ temp<-link(theta=beta, j=y[i]+k)$h #temp<-link(theta1=beta[1],theta2=beta[2],theta3=beta[3],theta4=beta[4],j=y[i]+k)$h si=si+(-1)^k*choose(myi, k)*temp } si=choose(m[i],y[i])*si loglik<- loglik+ log(ifelse(siF1[m]) Y[i]= m else{ for (k in 2:m) { if (U[i]>F1[k-1] && U[i]<= F1[k]) Y[i]=k-1 } } } return(Y) }