Program g02jafe

!     G02JAF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

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