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 ..
Continue
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 ..
Continue
! 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 27.3 Release. NAG Copyright 2021.
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