    PROGRAM g02jafe

!      G02JAF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g02jaf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: reml, tol
       INTEGER                         :: cwid, df, fint, i, ifail, j, k, l,   &
                                          lb, lddat, maxit, n, ncol, nff, nfv, &
                                          nrf, nrv, nvpr, rint, svid, warn, yvid
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: b(:), dat(:,:), gamma(:), se(:)
       INTEGER, ALLOCATABLE            :: fvid(:), levels(:), rvid(:), vpr(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          max
!      .. Executable Statements ..
       WRITE (nout,*) 'G02JAF Example Program Results'
       WRITE (nout,*)

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

!      Read in the problem size
       READ (nin,*) n, ncol, nfv, nrv, nvpr

       ALLOCATE (levels(ncol),fvid(nfv),rvid(nrv))

!      Read in number of levels for each variable
       READ (nin,*) levels(1:ncol)

!      Read in model information
       READ (nin,*) yvid, fvid(1:nfv), rvid(1:nrv), svid, cwid, fint, rint

!      If no subject specified, then ignore RINT
       IF (svid==0) THEN
          rint = 0
       END IF

!      Calculate LB
       lb = rint
       DO i = 1, nrv
          lb = lb + levels(rvid(i))
       END DO
       IF (svid/=0) THEN
          lb = lb*levels(svid)
       END IF
       lb = lb + fint
       DO i = 1, nfv
          lb = lb + max(levels(fvid(i))-1,1)
       END DO

       lddat = n
       ALLOCATE (vpr(nrv),dat(lddat,ncol),gamma(nvpr+2),b(lb),se(lb))

!      Read in the variance component flag
       READ (nin,*) vpr(1:nrv)

!      Read in the Data matrix
       READ (nin,*) (dat(i,1:ncol),i=1,n)

!      Read in the initial values for GAMMA
       READ (nin,*) gamma(1:(nvpr+rint))

!      Read in the maximum number of iterations
       READ (nin,*) maxit

!      Use default value for tolerance
       tol = 0.0E0_nag_wp

!      Fit the linear mixed effects regresion model
       ifail = 0
       CALL g02jaf(n,ncol,lddat,dat,levels,yvid,cwid,nfv,fvid,fint,nrv,rvid, &
          nvpr,vpr,rint,svid,gamma,nff,nrf,df,reml,lb,b,se,maxit,tol,warn, &
          ifail)

!      Display results
       IF (warn/=0) THEN
          WRITE (nout,*) 'Warning: At least one variance component was ', &
             'estimated to be negative and then reset to zero'
          WRITE (nout,*)
       END IF
       WRITE (nout,*) 'Fixed effects (Estimate and Standard Deviation)'
       WRITE (nout,*)
       k = 1
       IF (fint==1) THEN
          WRITE (nout,99999) 'Intercept             ', b(k), se(k)
          k = k + 1
       END IF
       DO i = 1, nfv
          DO j = 1, levels(fvid(i))
             IF (levels(fvid(i))==1 .OR. j/=1) THEN
                WRITE (nout,99995) 'Variable', i, ' Level', j, b(k), se(k)
                k = k + 1
             END IF
          END DO
       END DO

       WRITE (nout,*)
       WRITE (nout,*) 'Random Effects (Estimate and Standard', ' Deviation)'
       WRITE (nout,*)
       IF (svid==0) THEN
          DO i = 1, nrv
             DO j = 1, levels(rvid(i))
                WRITE (nout,99995) 'Variable', i, ' Level', j, b(k), se(k)
                k = k + 1
             END DO
          END DO
       ELSE
          DO l = 1, levels(svid)
             IF (rint==1) THEN
                WRITE (nout,99998) 'Intercept for Subject Level', l, &
                   '         ', b(k), se(k)
                k = k + 1
             END IF
             DO i = 1, nrv
                DO j = 1, levels(rvid(i))
                   WRITE (nout,99997) 'Subject Level', l, ' Variable', i, &
                      ' Level', j, b(k), se(k)
                   k = k + 1
                END DO
             END DO
          END DO
       END IF

       WRITE (nout,*)
       WRITE (nout,*) ' Variance Components'
       WRITE (nout,99996) (i,gamma(i),i=1,nvpr+rint)

       WRITE (nout,*)
       WRITE (nout,99994) 'SIGMA^2     = ', gamma(nvpr+rint+1)
       WRITE (nout,99994) '-2LOG LIKE  = ', reml
       WRITE (nout,99993) 'DF          = ', df

99999  FORMAT (1X,A,2F10.4)
99998  FORMAT (1X,A,I4,A,2F10.4)
99997  FORMAT (1X,3(A,I4),2F10.4)
99996  FORMAT (1X,I4,F10.4)
99995  FORMAT (1X,2(A,I4),2F10.4)
99994  FORMAT (1X,A,F10.4)
99993  FORMAT (1X,A,I16)
    END PROGRAM g02jafe
