! Last change: DOS 30 Jul 2000 2:04 pm ! *** copyright 2000 *** ! *** filename chex99.f95 *** John F. Monahan ** ! ********************** program chex99 ! Example 9.9 -- Poisson Regression using Frome's data on lung ! cancer death rates by age group and smoking rate group ! See E.L. Frome (1983) Biometrics, pp.665-674 ! implicit none integer, parameter :: ni = 9 integer, parameter :: nj = 7 integer, parameter :: np = 15 integer, dimension(ni,nj) :: nij,yij integer i,j,ii,jj,k,nn,nd,krep,idet real, dimension(np) :: beta,dell real, dimension( (np*(np+1))/2 ) :: del2l real lam,d,xb,rll,det ! ! i = 1,...,9=ni age group ! j = 1,...,7=nj smoking rate group where nonsmokers are j=1 ! ! nij(i,j) number in each i,j category ! yij(i,j) number of deaths in category ! 20 format(2i2,i6,i3) 21 format(' number of patients',i8,' deaths',i8) 22 format(8f10.5) 24 format(/' iteration ',i3,' log-likelihood',f14.7) 25 format(' new estimates, 9 age group effects, 6 smoke rates') 26 format(' update step ') 27 format(' gradient -- dell ') ! ! read in the data ! open( unit=5, file='frome.dat') open( unit=6, file='chex99.out') ! nn = 0 nd = 0 do j = 1,nj do i = 1,ni read(5,20) ii,jj,nij(i,j),yij(i,j) nn = nn + nij(i,j) nd = nd + yij(i,j) end do ! loop on i end do ! loop on j write(6,21) nn,nd ! set initial rate values lam = real(nd)/real(nn) lam = log(lam) beta = 0. ! initialize most to zero do i = 1,ni beta(i) = lam end do ! loop on i write(6,25) write(6,22) beta ! start looping do krep = 1,12 ! ! create dell and del2l ! ! initialize to zero rll = 0. dell = 0. del2l = 0. ! go through data do j = 1,nj do i = 1,ni xb = beta(i) ! set beta(i) = alpha(i) age group effects if( j .ne. 1 ) xb = xb + beta(ni-1+j) ! set beta(ni-1+j) = delta(j) smoke rate group ! delta(1)= 0 nonsmokers if( xb .gt. -120. ) then lam = exp(xb) else lam = 0. end if ! ( xb .gt. -120. ) ! first log likelihood ! Poisson rate is nij*lam rll = rll + yij(i,j)*log(nij(i,j)*lam) - nij(i,j)*lam ! d = yij(i,j) - nij(i,j)*lam ! dell dell(i) = dell(i) + d if( j .ne. 1 ) dell(ni-1+j) = dell(ni-1+j) + d ! del2l in symmetric form k = (i*(i+1))/2 del2l(k) = del2l(k) + nij(i,j)*lam if( j .gt. 1 ) then k = ((ni-2+j)*(ni-1+j))/2 + i del2l(k) = del2l(k) + nij(i,j)*lam k = ((ni-1+j)*(ni+j))/2 del2l(k) = del2l(k) + nij(i,j)*lam end if ! ( j .gt. 1 ) ! updated dell and del2l for this obs end do ! loop on i end do ! loop on j ! write rll and dell write(6,24) krep,rll write(6,27) write(6,22) dell ! call chlzky(del2l,np,det,idet) ! solve for Newton step call chlzhi(del2l,np,dell) call chlzih(del2l,np,dell) ! write out next step write(6,26) write(6,22) dell ! move beta - notice sign adjustments do k = 1,np beta(k) = beta(k) + dell(k) end do ! loop on k write(6,25) write(6,22) beta ! end of loop end do ! loop on krep stop end program chex99 SUBROUTINE ADJUST(D,I) ! NORMALIZES D WHILE KEEPING CONSTANT VALUE OF D * (2**I) INTEGER, PARAMETER :: IBIG = -2147483644 REAL, INTENT (IN OUT) :: D INTEGER, INTENT(IN OUT) :: I IF ( D .GT. 0.0 ) THEN DO WHILE ( D .LT. 1.0 ) D = D*16. I = I-4 ENDDO DO WHILE ( D .GT. 16.0 ) D = D/16. I = I+4 ENDDO ELSE I = IBIG ! IF D < 0 THEN I = - BIG ENDIF RETURN END SUBROUTINE ADJUST SUBROUTINE CHLZKY(A,N,DET,IDET) ! CHLSKY COMPUTES THE CHOLESKY (SQUARE-ROOT) FACTORIZATION ! A = L * ( L - TRANSPOSE ) WHERE L IS LOWER TRIANGULAR ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! A IS OVERWRITTEN BY L IN ITS LOWER TRIANGULAR PART ! SUBROUTINE ADJUST KEEPS DETERMINANT FROM EXPLODING USING ! 1 LE DET LE 16 AND DETERMINANT OF A IS DET*2**IDET ! ! J F MONAHAN (DEC,1983) DEPT OF STAT, N C S U, RALEIGH, N C 27650 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( (N*(N+1))/2 ), INTENT(IN OUT) :: A REAL, INTENT(OUT) :: DET INTEGER, INTENT(OUT) :: IDET REAL T,S INTEGER I,J,K,IM1,JM1 ! INITIALIZE FOR DETERMINANT DET = 1. IDET = 0 ! DO I-TH ROW DO I = 1,N T = A( (I*(I+1))/2 ) ! FIRST ROW IS TRIVIAL IF(I .GT. 1) THEN IM1 = I-1 DO J=1,IM1 ! WORK ON (I,J)-TH ELEMENT S = A( (I*(I-1))/2 + J ) IF(J .GT. 1) THEN JM1 = J-1 DO K = 1,JM1 S = S - A( (J*(J-1))/2 + K ) * A( (I*(I-1))/2 + K ) ENDDO ! LOOP ON K ENDIF A( (I*(I-1))/2 + J )= S / A( (J*(J+1))/2 ) ! WORK ON DIAGONAL ELEMENT T = T - A( (I*(I-1))/2 + J ) * A( (I*(I-1))/2 + J ) ENDDO ! LOOP ON J ! FINISHED WITH LOOP ON J ! UPDATE DETERMINANT WITH DIAGONAL ENDIF DET = DET * T ! KEEP DET FROM EXPLODING CALL ADJUST(DET,IDET) ! ABORT IF DET IS NEGATIVE OR ZERO IF(IDET.LT.-2000000000) RETURN ! FINISH A POSITIVE DIAGONAL A( (I*(I+1))/2 ) = SQRT(T) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLZKY SUBROUTINE CHLZHI(A,N,Y) ! ! OVERWRITES Y WITH INVERSE( L ) * Y ! WHERE L IS A LOWER TRIANGULAR MATRIX (N BY N) STORED IN A ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! TO COMPUTE QUADRATIC FORM IN THE INVERSE OF A POSITIVE DEFINITE ! MATRIX A, THAT IS Y' INVERSE( A ) Y THEN ! FIRST: FACTOR A BY CHLSKY(A,N,DET,IDET) ! THIS COMPUTES LOWER TRIANGULAR L SO THAT A = L * L' ! SECOND: COMPUTE INVERSE(L) * Y BY CHLSHI(A,N,Y) ! THIRD: COMPUTE THE SUM OF SQUARES OF THE ELEMENTS OF Y ! ! J F MONAHAN (JULY,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( (N*(N+1))/2 ), INTENT(IN) :: A REAL, DIMENSION( N ), INTENT(IN OUT) :: Y INTEGER I,J,IM1 ! ! SINCE THE MATRIX IS LOWER TRIANGULAR, SOLVE FROM THE TOP DOWN ! DO I = 1,N ! BRANCH AROUND ON THE FIRST ONE IF( I .GT. 1 ) THEN IM1 = I - 1 DO J = 1,IM1 Y(I) = Y(I) - A( (I*(I-1))/2 + J )*Y(J) ENDDO ENDIF Y(I) = Y(I) / A( (I*(I+1))/2 ) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLZHI SUBROUTINE CHLZIH(A,N,Y) ! *** NOTICE TRANSPOSE *** ! OVERWRITES Y WITH INVERSE( L' ) * Y ! WHERE L IS A LOWER TRIANGULAR MATRIX (N BY N) STORED IN A ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! TO COMPUTE THE PRODUCT OF THE INVERSE OF A POSITIVE DEFINITE MATRIX ! AND A VECTOR Y, THAT IS INVERSE( A ) * Y THEN ! FIRST: FACTOR A BY CHLSKY(A,N,DET,IDET) ! THIS COMPUTES LOWER TRIANGULAR L SO THAT A = L * L' ! SECOND: COMPUTE INVERSE(L) * Y BY CHLSHI(A,N,Y) ! THIRD: COMPUTE INVERSE(L') * ( INVERSE(L)*Y ) BY CHLSIH(A,N,Y) ! ! J F MONAHAN (JULY,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( (N*(N+1))/2 ), INTENT(IN) :: A REAL, DIMENSION( N ), INTENT(IN OUT) :: Y INTEGER I,II,IP1,J ! ! SINCE THE MATRIX IS UPPER TRIANGULAR, SOLVE FROM THE BOTTOM UP ! ! DO THE LAST ONE SEPARATELY Y(N) = Y(N) / A( (N*(N+1))/2 ) ! IF N = 1 THEN QUIT ELSE START LOOPING IF( N .EQ. 1 ) RETURN DO II = 2,N ! I BELOW NOW GOES FROM N-1 DOWN TO 1 I = N + 1 - II IP1 = I + 1 DO J = IP1,N Y(I) = Y(I) - A( (J*(J-1))/2 + I )*Y(J) ENDDO ! LOOP ON J Y(I) = Y(I) / A( (I*(I+1))/2 ) ENDDO ! LOOP ON II,I RETURN END SUBROUTINE CHLZIH ! *** end of filename chex99.f95 *********************