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

NAG FL Interface Introduction
Example description
!   G02JEF Example Program Text
!   Mark 30.0 Release. NAG Copyright 2024.

    Module g02jefe_mod

!     G02JEF Example Program Module:
!            Parameters and Utility Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: print_results
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6
    Contains
      Subroutine print_results(n,nff,nlsv,nrf,fixed,lfixed,nrndm,rndm,ldrndm,  &
        nvpr,vpr,lvpr,gamma,effn,rnkx,ncov,lnlike,lb,id,ldid,b,se)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: lnlike
        Integer, Intent (In)           :: effn, lb, ldid, ldrndm, lfixed,      &
                                          lvpr, n, ncov, nff, nlsv, nrf,       &
                                          nrndm, nvpr, rnkx
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: b(lb), gamma(nvpr+1), se(lb)
        Integer, Intent (In)           :: fixed(lfixed), id(ldid,lb),          &
                                          rndm(ldrndm,nrndm), vpr(lvpr)
!       .. Local Scalars ..
        Integer                        :: aid, i, k, l, ns, nv, p, pb, tb,     &
                                          tdid, vid
        Character (120)                :: pfmt, tfmt
!       .. Executable Statements ..
!       Display the output
        Write (nout,*) 'Number of observations (N)                    = ', n
        Write (nout,*) 'Number of random factors (NRF)                = ', nrf
        Write (nout,*) 'Number of fixed factors (NFF)                 = ', nff
        Write (nout,*) 'Number of subject levels (NLSV)               = ',     &
          nlsv
        Write (nout,*) 'Rank of X (RNKX)                              = ',     &
          rnkx
        Write (nout,*) 'Effective N (EFFN)                            = ',     &
          effn
        Write (nout,*) 'Number of nonzero variance components (NCOV) = ', ncov

        Write (nout,99990) 'Parameter Estimates'
        tdid = nff + nrf*nlsv

        If (nrf>0) Then
          Write (nout,*)
          Write (nout,99990) 'Random Effects'
        End If
        pb = -999
        pfmt = ' '
        Do k = 1, nrf*nlsv
          tb = id(1,k)
          If (tb/=-999) Then
            vid = id(2,k)
            nv = rndm(1,tb)
            ns = rndm(3+nv,tb)
            Write (tfmt,*)(id(3+l,k),l=1,ns)
            If (pb/=tb .Or. tfmt/=pfmt) Then
              If (k/=1) Then
                Write (nout,*)
              End If
              Write (nout,99991) ' Subject: ', ('Variable ',rndm(3+nv+l,tb),   &
                ' (Level ',id(3+l,k),')',l=1,ns)
            End If
            If (vid==0) Then
!             Intercept
              Write (nout,99994) b(k), se(k)
            Else
!             VID'th variable specified in RNDM
              aid = rndm(2+vid,tb)
              If (id(3,k)==0) Then
                Write (nout,99992) aid, b(k), se(k)
              Else
                Write (nout,99993) aid, id(3,k), b(k), se(k)
              End If
            End If
            pfmt = tfmt
            pb = tb
          End If
        End Do

        If (nff>0) Then
          Write (nout,*)
          Write (nout,99990) 'Fixed Effects'
        End If
        Do k = nrf*nlsv + 1, tdid
          If (vid/=-999) Then
            vid = id(2,k)
            If (vid==0) Then
!             Intercept
              Write (nout,99997) b(k), se(k)
            Else
!             VID'th variable specified in FIXED
              aid = fixed(2+vid)
              If (id(3,k)==0) Then
                Write (nout,99995) aid, b(k), se(k)
              Else
                Write (nout,99996) aid, id(3,k), b(k), se(k)
              End If
            End If
          End If
        End Do

        Write (nout,*)
        Write (nout,*) 'Variance Components'
        Write (nout,*) ' Estimate      Parameter        Subject'
        Do k = 1, nvpr
          Write (nout,99999,Advance='NO') gamma(k)
          p = 0
          Do tb = 1, nrndm
            nv = rndm(1,tb)
            ns = rndm(3+nv,tb)
            If (rndm(2,tb)==1) Then
              p = p + 1
              If (vpr(p)==k) Then
                Write (nout,99988,Advance='NO')(rndm(3+nv+l,tb),l=1,ns)
              End If
            End If
            Do i = 1, nv
              p = p + 1
              If (vpr(p)==k) Then
                Write (nout,99989,Advance='NO') rndm(2+i,tb),                  &
                  (rndm(3+nv+l,tb),l=1,ns)
              End If
            End Do
          End Do
          Write (nout,*)
        End Do
        Write (nout,*)
        Write (nout,99998) 'SIGMA**2         = ', gamma(nvpr+1)
        Write (nout,99998) '-2LOG LIKELIHOOD = ', lnlike

        Return
99999   Format (1X,F10.5,5X)
99998   Format (1X,A,F15.5)
99997   Format (3X,'Intercept',20X,F10.4,1X,F10.4)
99996   Format (3X,'Variable ',I2,' (Level ',I2,')',7X,F10.4,1X,F10.4)
99995   Format (3X,'Variable ',I2,18X,F10.4,1X,F10.4)
99994   Format (5X,'Intercept',18X,F10.4,1X,F10.4)
99993   Format (5X,'Variable ',I2,' (Level ',I2,')',5X,F10.4,1X,F10.4)
99992   Format (5X,'Variable ',I2,16X,F10.4,1X,F10.4)
99991   Format (1X,A,4(A,I2,A,I2,A,1X))
99990   Format (1X,A)
99989   Format (1X,'Variable',1X,I2,5X,'Variables',1X,100(I2,1X))
99988   Format (1X,'Intercept',7X,'Variables',1X,100(I2,1X))
      End Subroutine print_results
    End Module g02jefe_mod
    Program g02jefe

!     G02JEF Example Main Program

!     .. Use Statements ..
      Use g02jefe_mod, Only: nin, nout, print_results
      Use nag_library, Only: g02jcf, g02jef, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: lnlike
      Integer                          :: effn, i, ifail, j, lb, ldcxx, ldcxz, &
                                          ldczz, lddat, ldid, ldrndm, lfixed,  &
                                          licomm, liopt, lrcomm, lropt, lvpr,  &
                                          lwt, n, ncol, ncov, nff, nl, nlsv,   &
                                          nrf, nrndm, nv, nvpr, nzz, rnkx
      Character (1)                    :: weight
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: b(:), cxx(:,:), cxz(:,:), czz(:,:),  &
                                          dat(:,:), gamma(:), rcomm(:),        &
                                          ropt(:), se(:), wt(:), y(:)
      Integer, Allocatable             :: fixed(:), icomm(:), id(:,:),         &
                                          iopt(:), levels(:), rndm(:,:),       &
                                          vpr(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..
      Write (nout,*) 'G02JEF Example Program Results'
      Write (nout,*)

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

!     Read in the problem size
      Read (nin,*) weight, n, ncol, nrndm, nvpr

!     Set LFIXED and LDRNDM to maximum value they could
!     be for this dataset
      lfixed = ncol + 2
      ldrndm = 3 + 2*ncol

      If (weight=='W' .Or. weight=='w') Then
        lwt = n
      Else
        lwt = 0
      End If
      lddat = n
      Allocate (dat(lddat,ncol),levels(ncol),y(n),wt(lwt),fixed(lfixed),       &
        rndm(ldrndm,nrndm))

!     Read in the number of levels associated with each of the
!     independent variables
      Read (nin,*) levels(1:ncol)

!     Read in the fixed part of the model
      Read (nin,*)

!     Number of variables
      Read (nin,*) fixed(1)
      nv = fixed(1)

!     Intercept
      Read (nin,*) fixed(2)

!     Variable IDs
      If (nv>0) Then
        Read (nin,*) fixed(3:(nv+2))
      End If

!     Read in the random part of the model
      lvpr = 0
      Do j = 1, nrndm
!       Skip header
        Read (nin,*)

!       Number of variables
        Read (nin,*) rndm(1,j)
        nv = rndm(1,j)

!       Intercept
        Read (nin,*) rndm(2,j)

!       Variable IDs
        If (nv>0) Then
          Read (nin,*)(rndm(i,j),i=3,nv+2)
        End If

!       Number of subject variables
        Read (nin,*) rndm(nv+3,j)
        nl = rndm(nv+3,j)

!       Subject variable IDs
        If (nl>0) Then
          Read (nin,*)(rndm(i,j),i=nv+4,nv+nl+3)
        End If
        lvpr = lvpr + rndm(2,j) + nv
      End Do

!     Read in the dependent and independent data
      If (lwt>0) Then
        Read (nin,*)(y(i),dat(i,1:ncol),wt(i),i=1,n)
      Else
        Read (nin,*)(y(i),dat(i,1:ncol),i=1,n)
      End If

      licomm = 2
      lrcomm = 0
      Allocate (icomm(licomm),rcomm(lrcomm))

!     Call the initialization routine once to get LRCOMM and LICOMM
      ifail = 0
      Call g02jcf(weight,n,ncol,dat,lddat,levels,y,wt,fixed,lfixed,nrndm,rndm, &
        ldrndm,nff,nlsv,nrf,rcomm,lrcomm,icomm,licomm,ifail)

!     Reallocate ICOMM and RCOMM
      licomm = icomm(1)
      lrcomm = icomm(2)
      Deallocate (icomm,rcomm)
      Allocate (icomm(licomm),rcomm(lrcomm))

!     Pre-process the data
      ifail = 0
      Call g02jcf(weight,n,ncol,dat,lddat,levels,y,wt,fixed,lfixed,nrndm,rndm, &
        ldrndm,nff,nlsv,nrf,rcomm,lrcomm,icomm,licomm,ifail)

!     Use the default options
      liopt = 0
      lropt = 0

!     Calculate LDID
      ldid = 0
      Do i = 1, nrndm
        nv = rndm(1,i)
        ldid = max(rndm(3+nv,i),ldid)
      End Do
      ldid = ldid + 3

      lb = nff + nrf*nlsv
      nzz = nrf*nlsv
      ldczz = nzz
      ldcxx = nff
      ldcxz = nff
      Allocate (vpr(lvpr),gamma(nvpr+1),id(ldid,lb),b(lb),se(lb),              &
        czz(ldczz,nzz),cxx(ldcxx,nff),cxz(ldcxz,nzz),iopt(liopt),ropt(lropt))

!     Read in VPR
      Read (nin,*) vpr(1:lvpr)

!     Read in GAMMA
      Read (nin,*) gamma(1:nvpr)

!     Perform the analysis
      ifail = -1
      Call g02jef(lvpr,vpr,nvpr,gamma,effn,rnkx,ncov,lnlike,lb,id,ldid,b,se,   &
        czz,ldczz,cxx,ldcxx,cxz,ldcxz,rcomm,icomm,iopt,liopt,ropt,lropt,ifail)
      If (ifail/=0) Then
        If (ifail<100) Then
          Go To 100
        End If
      End If

!     Display results
      Call print_results(n,nff,nlsv,nrf,fixed,lfixed,nrndm,rndm,ldrndm,nvpr,   &
        vpr,lvpr,gamma,effn,rnkx,ncov,lnlike,lb,id,ldid,b,se)

100   Continue

    End Program g02jefe