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

NAG FL Interface Introduction
Example description
    Module g02jhfe_mod
!     G02JHF Example Program Module
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: custom_labels, 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
      Subroutine custom_labels(what,hlmm,vnames,plab)
!       Get custom labels for the fixed or random parameter estimates
!       or estimates of the variance components as obtained from fitting
!       a linear mixed effects regression model using G02JHF

!       WHAT  - what to get labels for. This is passed to G22YDF, valid
!               values are 'X' for the fixed parameter estimates, 'Z' for
!               the random parameter estimates and 'V' for the variance
!               components
!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
        Use nag_library, Only: g22ydf
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: hlmm
        Character (*), Intent (In)     :: what
!       .. Array Arguments ..
        Character (*), Allocatable, Intent (Out) :: plab(:)
        Character (*), Intent (In)     :: vnames(:)
!       .. Local Scalars ..
        Integer                        :: i, ifail, ip, j, k, li, lisx, lplab, &
                                          lvinfo, nv, vi
        Character (1)                  :: intcpt
        Character (100)                :: line, term
!       .. Local Arrays ..
        Integer, Allocatable           :: isx(:), vinfo(:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: adjustl, trim
!       .. Executable Statements ..

!       Get the require array sizes
        lisx = 0
        lvinfo = 3
        lplab = 0
        Allocate (plab(lplab),vinfo(lvinfo),isx(lisx))
        ifail = 1
        Call g22ydf(what,c_null_ptr,hlmm,intcpt,ip,lisx,isx,lplab,plab,lvinfo, &
          vinfo,ifail)
        If (ifail/=102) Then
!         Trigger the full error
          ifail = 0
          Call g22ydf(what,c_null_ptr,hlmm,intcpt,ip,lisx,isx,lplab,plab,      &
            lvinfo,vinfo,ifail)
        End If

!       Reallocate the output arrays
        lplab = vinfo(2)
        lvinfo = vinfo(3)
        Deallocate (plab,vinfo)
        Allocate (plab(lplab),vinfo(lvinfo))

!       Generate the labelling information
        Call g22ydf(what,c_null_ptr,hlmm,intcpt,ip,lisx,isx,lplab,plab,lvinfo, &
          vinfo,ifail)

!       Because we call g22ydf with lplab > 0, plab will contain the default
!       labels created by g22ydf, so we could return at this point, however
!       we will use the information in vinfo to create our own custom labels

        k = 1
j_lp:   Do j = 1, ip
          nv = vinfo(k)
          k = k + 1

          If (nv<0) Then
!           estimate is just there as padding
            Cycle j_lp
          End If

          If (nv==0) Then
            line = 'Intercept'
          Else
            line = ''
            Do i = 1, nv
              vi = vinfo(k)
              li = vinfo(k+1)
              If (li/=0) Then
                Write (term,99999) trim(adjustl(vnames(vi))), li
              Else
                Write (term,99998) trim(adjustl(vnames(vi)))
              End If

              If (i==1) Then
                line = trim(line) // trim(term)
              Else
                line = trim(line) // '.' // trim(term)
              End If

!             we are ignoring the contrast identifier when constructing
!             these labels
              k = k + 3
            End Do
          End If

          plab(j) = line
        End Do j_lp

        Return
99999   Format (A,' (Lvl ',I0,')')
99998   Format (A)
      End Subroutine custom_labels
    End Module g02jhfe_mod

!   G02JHF Example Program Text
!   Mark 28.5 Release. NAG Copyright 2022.
    Program g02jhfe
!     .. Use Statements ..
      Use g02jhfe_mod, Only: custom_labels, nin, nout, read_line
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nag_library, Only: g02jff, g02jhf, g22yaf, g22ybf, g22zaf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Type (c_ptr)                     :: hddesc, hfixed, hlmm
      Real (Kind=nag_wp)               :: lnlike
      Integer                          :: effn, fnlsv, i, ifail, 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)                    :: 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(:)
      Character (200), Allocatable     :: vplab(:), xplab(:), zplab(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: adjustl, size, trim
!     .. Executable Statements ..
      Write (nout,*) 'G02JHF Example Program Results'
      Write (nout,*)

      hfixed = c_null_ptr
      hddesc = c_null_ptr
      hlmm = 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

      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),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 custom_labels('Z',hlmm,vnames,zplab)

        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 random parameter estimates
        Call custom_labels('X',hlmm,vnames,xplab)

        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 custom_labels('V',hlmm,vnames,vplab)

        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,A47,1X,F10.4,5X,F10.4)
99996 Format (F10.4)
99995 Format (A10,40X,A10,5X,A)
99994 Format (A10,40X,A10)
    End Program g02jhfe