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

NAG FL Interface Introduction
Example description
!   G22YDF Example Program Text
!   Mark 28.5 Release. NAG Copyright 2022.

    Module g22ydfe_mod
!     G22YDF Example Program Module:
!            Parameters and User-defined Routines

!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: fit_lm, 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 fit_lm(hform,intcpt,nobs,mx,x,isx,ip,y,plab)
!       Perform a multiple linear regression using G02DAF

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
        Use nag_library, Only: g02daf, g22znf, nag_wp
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: hform
        Integer, Intent (In)           :: ip, mx, nobs
        Character (*), Intent (In)     :: intcpt
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x(:,:), y(:)
        Integer, Intent (In)           :: isx(:)
        Character (*), Intent (In)     :: plab(:)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: rss, rvalue, tol
        Integer                        :: i, idf, ifail, irank, ivalue, ldq,   &
                                          ldx, lwt, optype
        Logical                        :: svd
        Character (200)                :: cvalue
        Character (1)                  :: weight
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: b(:), cov(:), h(:), p(:), q(:,:),   &
                                          res(:), se(:), wk(:), wt(:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: repeat, size, trim
!       .. Executable Statements ..

        ldx = size(x,1)

!       We are assuming un-weighted data
        weight = 'U'
        lwt = 0
        ldq = nobs
        Allocate (b(ip),cov((ip*ip+ip)/2),h(nobs),p(ip*(ip+                    &
          2)),q(ldq,ip+1),res(nobs),se(ip),wk(ip*ip+5*(ip-1)),wt(lwt))

!       Use suggested value for tolerance
        tol = 0.000001E0_nag_wp

!       Fit a regression model
        ifail = 0
        Call g02daf(intcpt,weight,nobs,x,ldx,mx,isx,ip,y,wt,rss,idf,b,se,cov,  &
          res,h,q,ldq,svd,irank,p,tol,wk,ifail)

!       Get the formula for the model being fit
        ifail = 0
        Call g22znf(hform,'Formula',ivalue,rvalue,cvalue,optype,ifail)

!       Display the results
        Write (nout,*) 'Model: ', trim(cvalue)
        Write (nout,*) '                               Parameter   Standard'
        Write (nout,*) 'Coefficients                   Estimate      Error'
        Write (nout,*) repeat('-',51)
        Do i = 1, ip
          Write (nout,99997) plab(i), b(i), se(i)
        End Do
        Write (nout,*) repeat('-',51)
        Write (nout,99998) 'Residual sum of squares = ', rss
        Write (nout,99999) 'Degrees of freedom      = ', idf

        Return
99999   Format (1X,A,I9)
99998   Format (1X,A,F9.4)
99997   Format (1X,A30,1X,F7.3,5X,F7.3)
      End Subroutine fit_lm
    End Module g22ydfe_mod

    Program g22ydfe

!     .. Use Statements ..
      Use g22ydfe_mod, Only: fit_lm, nin, nout, read_line
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nag_library, Only: g22yaf, g22ybf, g22ycf, g22ydf, g22zaf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Type (c_ptr)                     :: hddesc, hform, hxdesc
      Integer                          :: i, ifail, ip, lddat, ldx, lisx,      &
                                          lplab, lvinfo, lvnames, mx, nobs,    &
                                          nvar, sddat, sdx
      Character (200)                  :: formula
      Character (1)                    :: intcpt, what
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: dat(:,:), x(:,:), y(:)
      Real (Kind=nag_wp)               :: tx(0,0)
      Integer, Allocatable             :: isx(:), levels(:), vinfo(:)
      Integer                          :: tisx(0), tvinfo(3)
      Character (50), Allocatable      :: plab(:), vnames(:)
      Character (1)                    :: tplab(0)
!     .. Executable Statements ..

      Write (nout,*) 'G22YDF Example Program Results'
      Write (nout,*)

      hform = c_null_ptr
      hddesc = c_null_ptr
      hxdesc = 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,*) nobs, 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,nobs,nvar,levels,lvnames,vnames,ifail)

!     Read in the data matrix and response variable
      lddat = nobs
      sddat = nvar
      Allocate (dat(lddat,sddat),y(nobs))
      Read (nin,*)(dat(i,1:nvar),y(i),i=1,nobs)

!     Read in the formula for the full model, remove comments and parse it
      Call read_line(formula)
      ifail = 0
      Call g22yaf(hform,formula,ifail)

!     Calculate the size of the design matrix
      ldx = 0
      sdx = 0
      ifail = 1
      Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,tx,ldx,sdx,mx,ifail)
      If (ifail/=81 .And. ifail/=82 .And. ifail/=91 .And. ifail/=92) Then
!       Redisplay any error messages not related to X being too small
        ifail = 0
        Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,tx,ldx,sdx,mx,ifail)
      End If

!     Generate the design matrix
      ldx = nobs
      sdx = mx
      Allocate (x(ldx,sdx))
      ifail = 0
      Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,x,ldx,sdx,mx,ifail)

      what = 'X'

!     Get size of output arrays
      lvinfo = 3
      lisx = 0
      lplab = 0
      ifail = 1
      Call g22ydf(what,hform,hxdesc,intcpt,ip,lisx,tisx,lplab,tplab,lvinfo,    &
        tvinfo,ifail)
      If (ifail/=102) Then
!       Redisplay the error message
        ifail = 0
        Call g22ydf(what,hform,hxdesc,intcpt,ip,lisx,tisx,lplab,tplab,lvinfo,  &
          tvinfo,ifail)
      End If

!     Allocate output arrays
      lisx = tvinfo(1)
      lplab = tvinfo(2)
!     We don't need VINFO as we are using labels in PLAB
      lvinfo = 0
      Allocate (isx(lisx),plab(lplab),vinfo(lvinfo))

!     Get the ISX flag and parameter labels
      ifail = 0
      Call g22ydf(what,hform,hxdesc,intcpt,ip,lisx,isx,lplab,plab,lvinfo,      &
        vinfo,ifail)

!     Fit the full model and print the results
      Call fit_lm(hform,intcpt,nobs,mx,x,isx,ip,y,plab)

!     Read in the formula for the sub-model, remove comments and parse it
      Call read_line(formula)
      ifail = 0
      Call g22yaf(hform,formula,ifail)

!     Get the ISX flag and parameter labels, as the new model has to be
!     a sub-model of the original one, the output arrays, ISX, PLAB and VINFO
!     can be reused as they will be of sufficient size
      ifail = 0
      what = 'S'
      Call g22ydf(what,hform,hxdesc,intcpt,ip,lisx,isx,lplab,plab,lvinfo,      &
        vinfo,ifail)

      Write (nout,*)
      Write (nout,*)

!     Fit the sub-model and print the results
      Call fit_lm(hform,intcpt,nobs,mx,x,isx,ip,y,plab)

!     Clean-up the G22 handles
      ifail = 0
      Call g22zaf(hform,ifail)
      Call g22zaf(hddesc,ifail)
      Call g22zaf(hxdesc,ifail)

      Deallocate (dat,x,y)
      Deallocate (isx,levels,vinfo)
      Deallocate (plab,vnames)

    End Program g22ydfe