Program g12bafe ! G12BAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. 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)) ! 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 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 = 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 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