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 29.3 Release. NAG Copyright 2023.
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