! G02JDF Example Program Text
! Mark 26.1 Release. NAG Copyright 2016.
Module g02jdfe_mod
! G02JDF Example Program Module:
! Parameters and User-defined 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
End If
pb = tb
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 g02jdfe_mod
Program g02jdfe
! G02JDF Example Main Program
! .. Use Statements ..
Use g02jdfe_mod, Only: nin, nout, print_results
Use nag_library, Only: g02jcf, g02jdf, 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,*) 'G02JDF 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 and intercept
Read (nin,*) rndm(1,j)
Read (nin,*) rndm(2,j)
nv = rndm(1,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 g02jdf(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 .And. ifail<100) Then
Go To 100
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 g02jdfe