PROGRAM g12bafe

!      G12BAF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g02gcf, g12baf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: dev, tol
       INTEGER                         :: i, idf, ifail, ip, ip1, iprint,      &
                                          irank, ldv, ldz, lisi, lomega, m,    &
                                          maxit, n, nd, ndmax, ns
       CHARACTER (1)                   :: offset
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: b(:), cov(:), omega(:), res(:),      &
                                          sc(:), se(:), sur(:,:), t(:), tp(:), &
                                          v(:,:), wk(:), y(:), z(:,:)
       INTEGER, ALLOCATABLE            :: ic(:), isi(:), isz(:), iwk(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          count, eoshift, log, max, real
!      .. Executable Statements ..
       WRITE (nout,*) 'G12BAF Example Program Results'
       WRITE (nout,*)

!      Skip heading in data file
       READ (nin,*)

!      Read in problem size
       READ (nin,*) n, m, ns, maxit, iprint, offset

       IF (offset=='Y' .OR. offset=='y') THEN
          lomega = n
       ELSE
          lomega = 0
       END IF
       IF (ns>0) THEN
          lisi = n
       ELSE
          lisi = 0
       END IF
       ldz = n
       ndmax = n
       ldv = n
       ALLOCATE (z(ldz,m),isz(m),t(n),ic(n),omega(lomega),isi(lisi),res(n), &
          tp(ndmax),sur(ndmax,max(ns,1)),iwk(2*n),y(n))

!      Read in the data
       IF (ns>0) THEN
          IF (lomega==0) THEN
             READ (nin,*) (t(i),z(i,1:m),ic(i),isi(i),i=1,n)
          ELSE
             READ (nin,*) (t(i),z(i,1:m),ic(i),isi(i),omega(i),i=1,n)
          END IF
       ELSE
          IF (lomega==0) THEN
             READ (nin,*) (t(i),z(i,1:m),ic(i),i=1,n)
          ELSE
             READ (nin,*) (t(i),z(i,1:m),ic(i),omega(i),i=1,n)
          END IF
       END IF

!      Read in the variable indication
       READ (nin,*) isz(1:m)

!      Calculate number of parameters in the model
       ip = count(isz(1:m)>0)

!      We are fitting a mean in the exponential model, so arrays used by G02GCF
!      need to be based on IP + 1
       ip1 = ip + 1
       ALLOCATE (b(ip1),se(ip1),sc(ip),cov(ip1*(ip1+1)/2),wk(ip1*(ip1+ &
          9)/2+n),v(ldv,ip1+7))

!      Specifiy tolerance to use
       tol = 5.0E-5_nag_wp

!      Get initial estimates by fitting an exponential model
       DO i = 1, n
          y(i) = 1.0E0_nag_wp - real(ic(i),kind=nag_wp)
          v(i,7) = log(t(i))
       END DO

!      Fit exponential model
       ifail = -1
       CALL g02gcf('L','M','Y','U',n,z,ldz,m,isz,ip1,y,res,0.0E0_nag_wp,dev, &
          idf,b,irank,se,cov,v,ldv,tol,maxit,0,0.0E0_nag_wp,wk,ifail)
       IF (ifail/=0) THEN
          IF (ifail<5) THEN
             GO TO 20
          END IF
       END IF

!      Check exponential model was of full rank
       IF (irank/=ip1) THEN
          WRITE (nout,*) ' WARNING: covariates not of full rank'
       END IF

!      Move all parameter estimates down one so as to drop the parameter
!      estimate for the mean.
       b = eoshift(b,1)

!      Fit Cox proportional hazards model
       ifail = -1
       CALL g12baf('No-offset',n,m,ns,z,ldz,isz,ip,t,ic,omega,isi,dev,b,se,sc, &
          cov,res,nd,tp,sur,ndmax,tol,maxit,iprint,wk,iwk,ifail)
       IF (ifail/=0) THEN
          IF (ifail<5) THEN
             GO TO 20
          END IF
       END IF

!      Display results
       WRITE (nout,*) ' Parameter      Estimate', '       Standard Error'
       WRITE (nout,*)
       WRITE (nout,99999) (i,b(i),se(i),i=1,ip)
       WRITE (nout,*)
       WRITE (nout,99998) ' Deviance = ', dev
       WRITE (nout,*)
       WRITE (nout,*) '    Time     Survivor Function'
       WRITE (nout,*)
       ns = max(ns,1)
       WRITE (nout,99997) (tp(i),sur(i,1:ns),i=1,nd)

20     CONTINUE

99999  FORMAT (I6,10X,F8.4,10X,F8.4)
99998  FORMAT (A,E13.4)
99997  FORMAT (F10.0,5X,F8.4)
    END PROGRAM g12bafe