! G02JEF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module g02jefe_mod ! G02JEF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. Use nag_library, Only: nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: 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 non-zero 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 g02jefe_mod Program g02jefe ! G02JEF Example Main Program ! .. Use Statements .. Use nag_library, Only: g02jcf, g02jef, nag_wp Use g02jefe_mod, Only: nin, nout, print_results ! .. 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 + 1 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 initialisation 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