Program g02jafe ! G02JAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: g02jaf, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: reml, tol Integer :: cwid, df, fint, i, ifail, j, k, l, & lb, lddat, maxit, n, ncol, nff, nfv, & nrf, nrv, nvpr, rint, svid, warn, yvid ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: b(:), dat(:,:), gamma(:), se(:) Integer, Allocatable :: fvid(:), levels(:), rvid(:), vpr(:) ! .. Intrinsic Procedures .. Intrinsic :: max ! .. Executable Statements .. Write (nout,*) 'G02JAF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) ! Read in the problem size Read (nin,*) n, ncol, nfv, nrv, nvpr Allocate (levels(ncol),fvid(nfv),rvid(nrv)) ! Read in number of levels for each variable Read (nin,*) levels(1:ncol) ! Read in model information Read (nin,*) yvid, fvid(1:nfv), rvid(1:nrv), svid, cwid, fint, rint ! If no subject specified, then ignore RINT If (svid==0) Then rint = 0 End If ! Calculate LB lb = rint Do i = 1, nrv lb = lb + levels(rvid(i)) End Do If (svid/=0) Then lb = lb*levels(svid) End If lb = lb + fint Do i = 1, nfv lb = lb + max(levels(fvid(i))-1,1) End Do lddat = n Allocate (vpr(nrv),dat(lddat,ncol),gamma(nvpr+2),b(lb),se(lb)) ! Read in the variance component flag Read (nin,*) vpr(1:nrv) ! Read in the Data matrix Read (nin,*)(dat(i,1:ncol),i=1,n) ! Read in the initial values for GAMMA Read (nin,*) gamma(1:(nvpr+rint)) ! Read in the maximum number of iterations Read (nin,*) maxit ! Use default value for tolerance tol = 0.0E0_nag_wp ! Fit the linear mixed effects regresion model ifail = 0 Call g02jaf(n,ncol,lddat,dat,levels,yvid,cwid,nfv,fvid,fint,nrv,rvid, & nvpr,vpr,rint,svid,gamma,nff,nrf,df,reml,lb,b,se,maxit,tol,warn,ifail) ! Display results If (warn/=0) Then Write (nout,*) 'Warning: At least one variance component was ', & 'estimated to be negative and then reset to zero' Write (nout,*) End If Write (nout,*) 'Fixed effects (Estimate and Standard Deviation)' Write (nout,*) k = 1 If (fint==1) Then Write (nout,99999) 'Intercept ', b(k), se(k) k = k + 1 End If Do i = 1, nfv Do j = 1, levels(fvid(i)) If (levels(fvid(i))==1 .Or. j/=1) Then Write (nout,99995) 'Variable', i, ' Level', j, b(k), se(k) k = k + 1 End If End Do End Do Write (nout,*) Write (nout,*) 'Random Effects (Estimate and Standard', ' Deviation)' Write (nout,*) If (svid==0) Then Do i = 1, nrv Do j = 1, levels(rvid(i)) Write (nout,99995) 'Variable', i, ' Level', j, b(k), se(k) k = k + 1 End Do End Do Else Do l = 1, levels(svid) If (rint==1) Then Write (nout,99998) 'Intercept for Subject Level', l, ' ', & b(k), se(k) k = k + 1 End If Do i = 1, nrv Do j = 1, levels(rvid(i)) Write (nout,99997) 'Subject Level', l, ' Variable', i, ' Level', & j, b(k), se(k) k = k + 1 End Do End Do End Do End If Write (nout,*) Write (nout,*) ' Variance Components' Write (nout,99996)(i,gamma(i),i=1,nvpr+rint) Write (nout,*) Write (nout,99994) 'SIGMA^2 = ', gamma(nvpr+rint+1) Write (nout,99994) '-2LOG LIKE = ', reml Write (nout,99993) 'DF = ', df 99999 Format (1X,A,2F10.4) 99998 Format (1X,A,I4,A,2F10.4) 99997 Format (1X,3(A,I4),2F10.4) 99996 Format (1X,I4,F10.4) 99995 Format (1X,2(A,I4),2F10.4) 99994 Format (1X,A,F10.4) 99993 Format (1X,A,I16) End Program g02jafe