NAG Library Manual, Mark 28.6
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
    Module g02jffe_mod
!     G02JFF Example Program Module
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: read_line
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6
    Contains
      Subroutine read_line(v1)
!       Read in a line from NIN and remove any comments

!       .. Scalar Arguments ..
        Character (*), Intent (Out)    :: v1
!       .. Local Scalars ..
        Integer                        :: pend
!       .. Intrinsic Procedures ..
        Intrinsic                      :: adjustl, index
!       .. Executable Statements ..

        Read (nin,'(A200)') v1
        pend = index(v1,'::')
        If (pend/=0) Then
          v1 = v1(1:pend-1)
        End If
        v1 = adjustl(v1)

        Return
      End Subroutine read_line
    End Module g02jffe_mod

!   G02JFF Example Program Text
!   Mark 28.6 Release. NAG Copyright 2022.
    Program g02jffe
!     .. Use Statements ..
      Use g02jffe_mod, Only: nin, nout, read_line
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nag_library, Only: g02jff, g02jhf, g22yaf, g22ybf, g22ydf, g22zaf,   &
                             g22zmf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Type (c_ptr)                     :: hddesc, hfixed, hlmm, hsform
      Real (Kind=nag_wp)               :: lnlike
      Integer                          :: effn, fnlsv, i, ifail, ip, lb,       &
                                          ldcxx, ldcxz, ldczz, lddat, licomm,  &
                                          lisx, lrcomm, lvinfo, lvnames, lwt,  &
                                          n, ncov, nff, nrf, nrndm, nvar,      &
                                          nvpr, nxx, nzz, rnkx, rnlsv, sddat
      Character (200)                  :: formula, line
      Character (1)                    :: intcpt, weight
!     .. Local Arrays ..
      Type (c_ptr), Allocatable        :: hrndm(:)
      Real (Kind=nag_wp), Allocatable  :: b(:), cxx(:,:), cxz(:,:), czz(:,:),  &
                                          dat(:,:), gamma(:), rcomm(:), se(:), &
                                          wt(:), y(:)
      Integer, Allocatable             :: icomm(:), isx(:), levels(:),         &
                                          vinfo(:)
      Character (50), Allocatable      :: vnames(:), vplab(:), xplab(:),       &
                                          zplab(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: adjustl, size, trim
!     .. Executable Statements ..
      Write (nout,*) 'G02JFF Example Program Results'
      Write (nout,*)

      hfixed = c_null_ptr
      hddesc = c_null_ptr
      hlmm = c_null_ptr
      hsform = c_null_ptr

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

!     Read in size of the data matrix and number of variable labels supplied
      Read (nin,*) weight, n, nvar, lvnames

!     Read in number of levels and names for the variables
      Allocate (levels(nvar),vnames(lvnames))
      Read (nin,*) levels(1:nvar)
      If (lvnames>0) Then
        Read (nin,*) vnames(1:lvnames)
      End If

!     Create a description of the data matrix
      ifail = 0
      Call g22ybf(hddesc,n,nvar,levels,lvnames,vnames,ifail)

!     Read in the data matrix and response variable
      lddat = n
      sddat = nvar
      If (weight=='W' .Or. weight=='w') Then
        lwt = n
      Else
        lwt = 0
      End If
      Allocate (dat(lddat,sddat),y(n),wt(lwt))
      If (lwt==0) Then
        Read (nin,*)(dat(i,1:nvar),y(i),i=1,n)
      Else
        Read (nin,*)(dat(i,1:nvar),y(i),wt(i),i=1,n)
      End If

!     Read in the formula for the fixed part of the model,
!     remove comments and parse it
      Call read_line(formula)
      ifail = 0
      Call g22yaf(hfixed,formula,ifail)

!     Read in the number of random statements
      Read (nin,*) nrndm

      Allocate (hrndm(nrndm))
      hrndm(:) = c_null_ptr

!     Read in the formula for the random parts of the model,
!     remove comments and parse them
      Do i = 1, nrndm
        Call read_line(formula)
        ifail = 0
        Call g22yaf(hrndm(i),formula,ifail)
      End Do

!     Get the size of the communication arrays
      licomm = 2
      lrcomm = 0
      Allocate (icomm(licomm),rcomm(lrcomm))
      ifail = 1
      Call g02jff(hlmm,hddesc,hfixed,nrndm,hrndm,weight,n,y,wt,dat,lddat,      &
        sddat,fnlsv,nff,rnlsv,nrf,nvpr,rcomm,lrcomm,icomm,licomm,ifail)
      If (ifail==191) Then
!       Allocate the communication arrays
        licomm = icomm(1)
        lrcomm = icomm(2)
        Deallocate (icomm,rcomm)
        Allocate (icomm(licomm),rcomm(lrcomm))
      End If

!     Pre-process the data
      ifail = -1
      Call g02jff(hlmm,hddesc,hfixed,nrndm,hrndm,weight,n,y,wt,dat,lddat,      &
        sddat,fnlsv,nff,rnlsv,nrf,nvpr,rcomm,lrcomm,icomm,licomm,ifail)
      If (ifail/=0 .And. ifail/=34 .And. ifail/=102) Then
!       terminate on an error
        Stop
      End If

!     Clear up the handles that are no longer required
      ifail = 0
      Call g22zaf(hddesc,ifail)
      Call g22zaf(hfixed,ifail)
      Do i = 1, nrndm
        Call g22zaf(hrndm(i),ifail)
      End Do

!     Specify that we want to fit the model by maximum likelihood (ML)
      ifail = 0
      Call g22zmf(hlmm,'Likelihood = ML',ifail)

      nzz = nrf*rnlsv
      nxx = nff*fnlsv
      lb = nxx + nzz

!     We don't want the output from CXX, CZZ and CXZ
      Allocate (czz(0,0),cxx(0,0),cxz(0,0))
      ldczz = size(czz,1)
      ldcxx = size(cxx,1)
      ldcxz = size(cxz,1)

!     Allocate the rest of the output arrays
      Allocate (gamma(nvpr+1),b(lb),se(lb))

!     Use MIVQUE for the initial values for gamma
      gamma(:) = -1.0_nag_wp

!     Fit the model
      ifail = -1
      Call g02jhf(hlmm,nvpr,gamma,effn,rnkx,ncov,lnlike,lb,b,se,czz,ldczz,cxx, &
        ldcxx,cxz,ldcxz,rcomm,icomm,ifail)
      If (ifail/=0 .And. ifail/=1001 .And. ifail/=1002 .And. ifail/=1003 .And. &
        ifail/=1004) Then
!       terminate if there is an error
        Stop
      End If

      lvinfo = 0
      lisx = 0
      Allocate (vinfo(lvinfo),xplab(nxx),zplab(nzz),vplab(nvpr),isx(lisx))

!     Print the results
      Write (nout,99998) 'Number of observations           = ', n
      Write (nout,99998) 'Total number of random factors   = ', nzz
      Write (nout,99998) 'Total number of fixed factors    = ', nxx
      Write (nout,99998) 'Rank of X                        = ', rnkx
      Write (nout,99998) 'Effective N                      = ', effn
      Write (line,99996) lnlike
      line = adjustl(line)
      Write (nout,99999) '-2Log Likelihood                 = ', trim(line)
      Write (line,99996) gamma(nvpr+1)
      line = adjustl(line)
      Write (nout,99999) 'Sigma**2                         = ', trim(line)

      Write (nout,99999) 'Parameter Estimates'
      If (nzz>0) Then
!       Get labels for the random parameter estimates
        Call g22ydf('Z',hsform,hlmm,intcpt,ip,lisx,isx,nzz,zplab,lvinfo,vinfo, &
          ifail)

        Write (nout,*)
        Write (nout,99999) 'Random Effects'
        Write (nout,99995) 'Parameter', 'Estimate', 'Standard Error'
        Do i = 1, nzz
          If (zplab(i)/='') Then
            Write (nout,99997) adjustl(zplab(i)), b(i), se(i)
          End If
        End Do
      End If

      If (nxx>0) Then
!       Get labels for the fixed parameter estimates
        Call g22ydf('X',hsform,hlmm,intcpt,ip,lisx,isx,nxx,xplab,lvinfo,vinfo, &
          ifail)

        Write (nout,*)
        Write (nout,99999) 'Fixed Effects'
        Write (nout,99995) 'Parameter', 'Estimate', 'Standard Error'
        Do i = 1, nxx
          Write (nout,99997) adjustl(xplab(i)), b(nzz+i), se(nzz+i)
        End Do
      End If

      If (nvpr>0) Then
!       Get labels for the variance component estimates
        Call g22ydf('V',hsform,hlmm,intcpt,ip,lisx,isx,nvpr,vplab,lvinfo,      &
          vinfo,ifail)

        Write (nout,*)
        Write (nout,99999) 'Variance Components'
        Write (nout,99994) 'Component', 'Estimate'
        Do i = 1, nvpr
          Write (*,99997) adjustl(vplab(i)), gamma(i)
        End Do
      End If

!     Clean-up the remaining G22 handles
      ifail = 0
      Call g22zaf(hlmm,ifail)

99999 Format (1X,A,A)
99998 Format (1X,A,I0)
99997 Format (1X,A25,2(5X,F10.4))
99996 Format (F10.4)
99995 Format (A10,22X,A10,5X,A)
99994 Format (A10,22X,A10)
    End Program g02jffe