Program g12bafe

!     G12BAF Example Program Text

!     Mark 27.2 Release. NAG Copyright 2021.

!     .. 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 Procedures ..
      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))

!     Specify 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 100
        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(1:ip1) = eoshift(b(1:ip1),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 100
        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)

100   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