#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
# eha.f
# intcovcb.h
# optioncb.h
# This archive created: Wed Jan 25 15:50:19 1995
export PATH; PATH=/bin:/usr/bin:$PATH
if test -f 'eha.f'
then
echo shar: "will not over-write existing file 'eha.f'"
else
cat << \SHAR_EOF > 'eha.f'
************************************************************************
* *
* eha (also known as wfsfit) Version 1.0 1-25-95 *
* *
* Program to generate replicate covariate data for each discrete *
* time unit (a single year in this implementation). The requisite *
* data from the separate World Fertility Study data files is copied *
* (for as many years as necessary) into one very big event history *
* covariate array: intcovar. *
* *
* Developer: Steven M. Lewis (slewis@stat.washington.edu). *
* This software is not formally maintained, but I will be happy to *
* hear from people who have problems with it, although I cannot *
* guarantee that all problems will be corrected. *
* *
* Permission is hereby granted to StatLib to redistribute this *
* software. The software can be freely used for non-commercial *
* purposes, and can be freely distributed for non-commercial *
* purposes only. The copyright is retained by the developer. *
* *
* Copyright 1995 Steven M. Lewis *
* *
* *
* References: *
* *
* Raftery, A.E., Lewis, S.M., Aghajanian, A. and Kahn, M.J. (1993). *
* "Event History Modeling of World Fertility Survey Data." Working *
* Paper No. 93-1, Center for Studies in Demography and Ecology, *
* University of Washington, Seattle, Wa. *
* This paper is available via the World Wide Web by linking to URL *
* http://www.stat.washington.edu/tech.reports/pub/tech.reports *
* and then selecting the "evenhist.ps" link. *
* This paper is also available via regular ftp using the following *
* commands: *
* ftp ftp.stat.washington.edu (or 128.95.17.34) *
* login as anonymous *
* enter your email address as your password *
* ftp> cd /pub/tech.reports *
* ftp> get evenhist.ps *
* ftp> quit *
* *
* Raftery, A.E., Lewis, S.M. and Aghajanian, A. (1994). "Demand or *
* Ideation? Evidence from the Iranian Marital Fertility Decline." *
* Working Paper No. 94-1, Center for Studies in Demography and *
* Ecology, University of Washington, Seattle, Wa. *
* This paper is available via the World Wide Web by linking to URL *
* http://www.stat.washington.edu/tech.reports/pub/tech.reports *
* and then selecting the "demidea.ps" link. *
* This paper is also available via regular ftp using the following *
* commands: *
* ftp ftp.stat.washington.edu (or 128.95.17.34) *
* login as anonymous *
* enter your email address as your password *
* ftp> cd /pub/tech.reports *
* ftp> get demidea.ps *
* ftp> quit *
* *
************************************************************************
program wfsfit
************************************************************************
* *
* The fortran commands defining the intcovar common block (intcovar *
* is the large internal event history covariates matrix) are in a *
* separate include file, 'intcovcb.h'. This include file also *
* contains definitions for four related parameters; the key one being *
* ppmaxcov which is set to the maximum number of covariates allocated *
* for the current implementation of these programs. Use of this *
* include file is not in accordance with the Fortran 77 Ansi *
* standard. *
* *
************************************************************************
include 'intcovcb.h'
************************************************************************
* *
* The fortran commands defining the option common block are in a *
* separate include file, 'optioncb.h'. Use of this include file is *
* not in accordance with the Fortran 77 Ansi standard. *
* *
************************************************************************
include 'optioncb.h'
************************************************************************
* *
* The user's parity selection is returned by the modelin subroutine *
* in the following two variables. Partype is a single character *
* which is either a blank or a '+'. Parint contains a number which *
* represents the first (chronologically earliest) parity to include *
* in the model. If partype is a blank, then only birth intervals *
* with the parity indicated by parint are to be included in the *
* model. If partype is a '+', then all birth intervals with parities *
* greater than or equal to parint are to be included in the model. *
* *
************************************************************************
character partype*1
integer parint
************************************************************************
* *
* Only a single record of the birth event history file will ever be *
* needed at any given point in this program. Therefore, no event *
* history vector need be declared. Instead, simply a scalar *
* character string, evenhist, is declared. It will be expedient to *
* use a read statement to substring the current birth event history *
* record into its separate subfields. The outputs from the read *
* statement are stored in the character strings located in the *
* eventcb common block. The data in this common block are also used *
* by the rowgen subroutine. Note that these declarations are listed *
* in the order in which the subfields exist in the evenhist *
* character string. *
* *
************************************************************************
common /eventcb/ intbegcm, intbegyr, intendcm, intendyr, closed,
+ parity, previous, prevlivm, prevlivy, breastfd, brestfln
character intbegcm*3
character intbegyr*2
character intendcm*3
character intendyr*2
character closed*1
character parity*2
character previous*2
character prevlivm*3
character prevlivy*2
character breastfd*2
character brestfln*2
************************************************************************
* *
* The covnames character vector will contain the names of all the *
* covariates to be included in this run of the model. It is also *
* used while writing the logistic regression results. The need to *
* check for an end-of-file explains why covnames contains an extra *
* "apparently unnecessary" element (see the comments in the modelin *
* subroutine for a more detailed explanation of this). *
* *
************************************************************************
character covnames(ppmaxcp1)*60
************************************************************************
* *
* The following variables hold various temporary values used in *
* different segments of this program (these are double precision or *
* character variables). *
* *
* fixed - vector of current estimated fixed effect parameters *
* random - vector of current estimated random effect parameters *
* trials - a vector containing the number of bernoulli trials *
* for each binomial sample (currently this is always *
* set to 1) *
* success - number of successes (births) during the current year *
* hessian - the hessian matrix for the likelihood *
* varmat - estimated fixed effect parameter covariance matrix *
* stddevs - standard deviations of all covariates in the model *
* means - sample means of all covariates in the model *
* llk - loglikelihood of the full model (ie, the model with *
* all the intcovar covariates included) *
* ll2 - -2 * the difference in the loglikelihood function *
* for the full model minus the model with only an *
* intercept term *
* ybar - average number of births per woman exposure-year *
* *
* province - vector holding province level data for all provinces *
* provuse - vector indicating which provinces are to be included *
* covterms - matrix containing parsed version of covariate names *
* covar - vector holding complete covariate data for all women *
* evenhist - character string containing all current event data *
* inpfile - global options input file name from command line *
* *
************************************************************************
double precision fixed(ppmaxcov) /ppmaxcov*0.0d0/
double precision random(4934) /4934*0.0d0/
double precision trials(93100)
double precision success(93100)
double precision hessian(ppmatmax)
double precision varmat(ppmaxcov,ppmaxcov)
double precision stddevs(ppmaxcov)
double precision means(ppmaxcov)
double precision llk
double precision ll2
double precision ybar
character province(0:23)*97
character provuse(0:22)*1
character covterms(ppmaxcov,3)*18
character covar(0:4934)*140
character evenhist*39
character inpfile*24
************************************************************************
* *
* Additional variables for storing various temporary values used in *
* different segments of this program (these are integer variables). *
* *
* womenids - women ids associated with each row of intcovar *
* womaxcnt - maximum number of random effect parameters available *
* provmax - maximum valid Iranian province number *
* maxiter - maximum number of iterations in the logit subroutine *
* maxorder - highest order interaction permitted in the model *
* termcnt - vector containing actual covariate interaction orders *
* argcount - total number of arguments passed in the command line *
* numcovar - number of columns in the input covariate matrix *
* duration - number of years since the woman's last birth *
* birthcnt - running total count of number of births *
* logitrc - error return code from the logit subroutine *
* intcnt - row index into intcovar for the current interval *
* iargc - Sun subroutine to get number of command line arguments *
* begyr - beginning of the current event (in Iranian years) *
* endyr - end of the current event (in Iranian years) *
* iyear - current Iranian calendar year for each exposure-year *
* yob - Iranian calendar year in which the woman was born *
* age - woman's current age for each exposure-year *
* *
************************************************************************
integer womenids(93100)
integer womaxcnt /4934/
integer provmax /22/
integer maxiter /16/
integer maxorder /3/
integer termcnt(ppmaxcov)
integer argcount
integer numcovar
integer duration
integer birthcnt
integer logitrc
integer intcnt
integer iargc
integer begyr
integer endyr
integer iyear
integer yob
integer age
************************************************************************
* *
* The following variables are intermediate work variables, do-loop *
* counters and similar such temporary subscripts and indices. *
* *
************************************************************************
integer ev
integer id
integer r1
integer c1
integer x1
integer rc
************************************************************************
* *
* Process whatever command line arguments were passed to wfsfit. *
* Currently at most one argument is expected. If present, the first *
* argument contains the name of the global options input file. If *
* not present, set the global options input file name to a ' '. *
* *
************************************************************************
argcount = iargc()
if (argcount.ge.1) then
call getarg(1,inpfile)
else
inpfile = ' '
end if
************************************************************************
* *
* Call optionin to initialize the options common block to all its *
* default values and, if the name of a global options input file was *
* passed as the first (and for now, only) command line argument, *
* replace the default values with any user specified options read *
* from the global options input file. *
* *
************************************************************************
call optionin(inpfile,rc)
if (rc.ne.0) then
write (0,*) 'optionin exited with a nonzero error code of', rc
go to 900
end if
************************************************************************
* *
* Enable my exception signal handler. However, since the rowgen and *
* logit subroutines generate a very large number of floating point *
* inexact exceptions, I had better disable just this exception type. *
* *
************************************************************************
if (ieee.eq.'T') call handlon('all-inexact',rc)
************************************************************************
* *
* Initialize the following two static common resident constants. *
* *
************************************************************************
maxint = 93100
maxcovar = ppmaxcov
************************************************************************
* *
* In order to initialize the parity selection variables, parint and *
* partype, and the vector of the names of the covariates to include *
* in this model run, call the modelin subroutine. This subroutine *
* copies the pertinent data from the model input file into these *
* variables. A zero return code indicates all went as expected. If *
* the return code was anything else, then modelin has already written *
* an explanatory error message and there is nothing left for us to do *
* except to terminate the program. *
* *
************************************************************************
call modelin(modelname,maxcovar,maxorder,covterms,covnames,
+ numcovar,termcnt,parint,partype,rc)
if (rc.ne.0) go to 800
************************************************************************
* *
* Call prvparse to convert the provset option, located in the options *
* common block, into the corresponding vector of 'T' or 'F' character *
* strings, provuse, for use by the chkevent subroutine. A zero *
* return code indicates all went as expected. If the return code was *
* anything else, then prvparse has already written an explanatory *
* error message and there is nothing left for us to do except to *
* terminate the program. *
* *
************************************************************************
call prvparse(provset,provmax,provuse,rc)
if (rc.ne.0) go to 800
************************************************************************
* *
* Call the provin subroutine to read the province level variables *
* into the province character string vector. Currently no errors are *
* returned via the return code, rc. *
* *
************************************************************************
call provin(provmax,province,rc)
************************************************************************
* *
* Call the covarin subroutine to read the separate woman covariate *
* file into the covariate character string vector. Currently no *
* errors are returned via the return code, rc. *
* *
************************************************************************
call covarin(womaxcnt,covar,rc)
************************************************************************
* *
* Now read the fertility event history file and obtain the correct *
* value of all database variables associated with the current birth *
* interval. Clearly not all of these are needed, but at least they *
* will all be available for easy reference. The chkevent subroutine *
* is called to do this. It initializes the provincb, covarcb and *
* eventcb common blocks. If the user's parity requirements are not *
* satisfied by this interval, or if any requisite variables are *
* missing or if there are invalid age, year or province variables *
* associated with this event, chkevent returns a nonzero error code *
* to us. If this is the case, this most recently read birth interval *
* should not be used while estimating the current model. *
* *
************************************************************************
open (unit=34,file='evenhist.output',status='old')
birthcnt = 0
intcnt = 0
do 500 ev=1,27500
read (34,'(a)',end=510) evenhist
call chkevent('complete',province,provmax,provuse,covar,
+ womaxcnt,evenhist,parint,partype,begyr,endyr,id,yob,rc)
if (rc.ne.0) go to 500
************************************************************************
* *
* Set years since last birth variable, duration, to 0. *
* *
************************************************************************
duration = 0
************************************************************************
* *
* Now add this birth interval to the arrays which will be passed to *
* the logit subroutine: intcovar, trials, and success. Add one *
* additional row to intcovar for each year contained in this *
* interval. In doing so, check that the number of cases has not *
* exceeded the maximum number of rows available in these 3 arrays. *
* *
************************************************************************
do 400 iyear=begyr,endyr
if (intcnt.ge.maxint) then
write (0,*) 'maximum number of cases has been exceeded'
close (unit=34)
go to 800
end if
intcnt = intcnt + 1
age = iyear - yob
************************************************************************
* *
* Call subroutine rowgen to calculate and store the entire set of *
* covariates for the current woman exposure-year. Only the row of *
* intcovar indexed by intcnt is altered by this subroutine. It uses *
* inputs from three common blocks, provincb, covarcb and eventcb, *
* along with its arguments to perform this task. If rowgen could not *
* complete its job, it returns a nonzero error return code via rc. *
* In this latter case, we had best just stop executing the program. *
* *
* Save the woman id number associated with the current exposure-year. *
* *
************************************************************************
call rowgen(covterms,maxcovar,termcnt,numcovar,maxint,intcnt,
+ duration,iyear,age,intcovar,rc)
if (rc.ne.0) then
close (unit=34)
go to 800
end if
womenids(intcnt) = id
************************************************************************
* *
* After the current row of intcovar is filled in, I must still set *
* the associated rows of trials and success. Each row represents a *
* single trial, so set the current row of trials to 1. Only the *
* final year of a closed interval actually represents a birth, so *
* while inside this loop assume that each year is not the final *
* year, and hence set the current row of success to 0. *
* *
* In case the current birth interval is not yet completed, increment *
* the year since last birth variable, duration, by 1. *
* *
************************************************************************
trials(intcnt) = 1.0d0
success(intcnt) = 0.0d0
duration = duration + 1
400 continue
************************************************************************
* *
* Set the success flag only for the year ending a closed interval. *
* Also update the running count of total number of births. *
* *
************************************************************************
if (closed.eq.'1') then
success(intcnt) = 1.0d0
birthcnt = birthcnt + 1
end if
500 continue
510 close (unit=34)
************************************************************************
* *
* In the case where the standardize option has been set to 'T'rue, *
* call the stdindep subroutine to apply a correlation transformation *
* to the intcovar matrix, but not to the first column since the logit *
* subroutine requires that the first column contain the intercept. *
* Note that this subroutine replaces the covariate values in the *
* intcovar matrix so the original values will no longer be directly *
* available. However the stdindep subroutine does return the sample *
* mean and standard deviation of each of the covariates so that the *
* original values may be derived from the transformed variables. *
* *
************************************************************************
if (standardize.eq.'T'.and.numcovar.gt.1) then
call stdindep(intcovar(1,2),maxint,(numcovar-1),intcnt,
+ intcovar(1,2),means(2),stddevs(2),rc)
if (rc.ne.0) then
write (0,*) 'stdindep exited with a nonzero error code of', rc
go to 800
end if
end if
************************************************************************
* *
* In the case where the intcdump option has been set to 'T'rue, call *
* the wrtintc subroutine to dump the entire contents of the intcovar *
* matrix to the dump file named in dumpname (which is located in the *
* options common block). *
* *
************************************************************************
if (intcdump.eq.'T') then
call wrtintc(intcovar,maxint,numcovar,intcnt,womenids,rc)
if (rc.ne.0) then
write (0,*) 'wrtintc exited with a nonzero error code of', rc
go to 800
end if
end if
************************************************************************
* *
* Initialize the parameter corresponding to the intercept to be *
* log(ybar/(1 - ybar)). *
* *
* Then, call the logit subroutine to do the calculations. *
* *
************************************************************************
ybar = dble(birthcnt) / dble(intcnt)
fixed(1) = dlog( ybar / (1.0d0 - ybar) )
call logit(fixed,ll2,hessian,success,trials,intcovar,maxint,
+ intcnt,numcovar,maxiter,llk,logitrc)
if (logitrc.ne.0) then
write (0,*) 'logit exited with a nonzero error code of',
+ logitrc
go to 800
end if
************************************************************************
* *
* Construct the entire estimated covariance matrix of the parameters. *
* *
************************************************************************
do 700 r1=1,numcovar
do 600 c1=r1,numcovar
call sssub(numcovar,r1,c1,x1,rc)
if (rc.ne.0) then
write (0,*) 'sssub exited with a nonzero error code of', rc
go to 800
end if
varmat(r1,c1) = hessian(x1)
varmat(c1,r1) = varmat(r1,c1)
600 continue
700 continue
************************************************************************
* *
* Write out the final results of the logistic regression. This is *
* now done by calling the wrtfixed subroutine. This subroutine will *
* also optionally write out the entire estimated correlation matrix *
* of the parameters. *
* *
************************************************************************
call wrtfixed(varmat,maxcovar,fixed,numcovar,covnames,parint,
+ partype,ll2,llk,intcnt,birthcnt,maxiter,logitrc,rc)
if (rc.ne.0) then
write (0,*) 'wrtfixed exited with a nonzero error code of', rc
go to 800
end if
************************************************************************
* *
* In the case where the residout option has been selected, call the *
* wrtresid subroutine to write the woman's id number, observed value, *
* fitted value and residual for each exposure-year to the file named *
* in residname (which is located in the options common block). *
* *
************************************************************************
if (residout.eq.'T') then
call wrtresid(intcovar,maxint,numcovar,intcnt,fixed,random,
+ womaxcnt,womenids,success,rc)
if (rc.ne.0) then
write (0,*) 'wrtresid exited with a nonzero error code of', rc
go to 800
end if
end if
************************************************************************
* *
* Finish up the program by cleaning up after ourselves. In *
* particular, disable my exception signal handler. *
* *
************************************************************************
800 if (ieee.eq.'T') call handloff('all',rc)
900 stop
end
************************************************************************
* *
* optionin 10-28-94 *
* *
* This subroutine reads the global options input file named in the *
* first argument and sets the variables in the options common block *
* according to what is read in from this file. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* filename = a character string containing the name of a file which *
* contains the global options to use during this run of *
* the Iran Fertility Survey modeling program. *
* *
* *
* Outputs: *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = unable to open the global options input file *
* named in the filename argument. *
* 8 = the oneparse subroutine returned with a nonzero *
* error return code. *
* 12 = the translat subroutine returned with a nonzero *
* error return code. *
* 16 = an unrecognizable option keyword was read from *
* the global options input file. *
* 20 = an invalid option value was encountered. For *
* example, something other than a 'T' or 'F' was *
* found where a 'T'rue/'F'alse indicator was *
* expected. *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine optionin(filename,r15)
character filename*(*)
integer r15
************************************************************************
* *
* The fortran commands defining the option common block are in a *
* separate include file, 'optioncb.h'. Use of this include file is *
* not in accordance with the Fortran 77 Ansi standard. *
* *
************************************************************************
include 'optioncb.h'
************************************************************************
* *
* The following variables hold various constants and temporary values *
* used in this subroutine such as loop counters and so forth. *
* *
************************************************************************
double precision realval
character ucase*26 /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
character lcase*26 /'abcdefghijklmnopqrstuvwxyz'/
character delimit*1 /'='/
character optterms(2)*60
character optrecrd*80
character keyword*18
character optval*60
character torf*1
integer termmax /2/
integer termcnt
integer intval
integer rc
************************************************************************
* *
* Initialize the variables in the options common block to their *
* default values. These may then be overridden by whatever is read *
* from the global options input file. *
* *
************************************************************************
optinname = filename
modelname = 'ifsmodel'
fixdname = 'fixdinit'
randname = 'randinit'
dumpname = 'intcdump.output'
logitname = 'logitresults'
gibbsname = 'gibbsresults'
residname = 'resids.output'
mlfrprioname = 'mlfrprio.output'
printopt = '6c'
provset = '0-22'
fixedmethod = 'metropolis'
fixedprior = 'normal'
fixedsupport = '-inf:inf'
randprior = 'normal'
varprior = 'inverted gamma'
fixedmargvar = dexp(1.0d0)
varshape = 0.47d0
varscale = 0.0226d0
firstvariance = 0.25d0
iteracnt = 75500
burnin = 500
priosize = 1000000
parmdps = 6
cseed = 43827
tseed = 1073
randtracecnt = 15
fixedcnt = 1
randcnt = 1
fixedndx = 2
womanid = 3911
ieee = 'T'
standardize = 'F'
intcdump = 'F'
residout = 'F'
gibbstrace = 'T'
gibbstime = 'T'
mlfrpriotrace = 'F'
fixedchk = 'F'
randchk = 'F'
varchk = 'F'
keep4707 = 'F'
longeropen = 'T'
************************************************************************
* *
* If the filename argument is ' ', the caller did not pass a global *
* options input file name in the invoking command line. In this case *
* we should use only the default values, so should return to the *
* calling routine immediately. *
* *
************************************************************************
if (filename.eq.' ') then
r15 = 0
return
end if
************************************************************************
* *
* Open the global options input file. If unable to do so, return to *
* the caller with an error code of 4. *
* *
************************************************************************
open (unit=20,file=filename,status='old',iostat=rc)
if (rc.ne.0) then
r15 = 4
return
end if
************************************************************************
* *
* Read the options from the global options input file. *
* *
************************************************************************
200 read (20,'(a)',end=700) optrecrd
************************************************************************
* *
* Parse each input option record into its keyword and associated *
* value. The keyword is separated from its value by an '=' sign, *
* which serves as the special character delimit variable when calling *
* the oneparse subroutine below. The oneparse subroutine also *
* left justifies (removes leading blanks) both the keyword and value. *
* If any blank lines are read from the global options input file, *
* these should be ignored. *
* *
* After parsing the input option record, translate the option keyword *
* into all lowercase letters. In this fashion the program will not *
* be sensitive to the case used to specify the option keyword. Note *
* that the option value is not translated so that the value's case is *
* significant (upper and lower case letters are interpreted as *
* different characters). *
* *
************************************************************************
call oneparse(optrecrd,delimit,termmax,optterms,termcnt,rc)
if (rc.ne.0) then
write (0,*) 'oneparse exited with a nonzero error code of', rc
r15 = 8
go to 800
end if
if (termcnt.lt.1) go to 200
call translat(ucase,lcase,optterms(1),1,keyword,rc)
if (rc.ne.0) then
write (0,*) 'translat exited with a nonzero error code of', rc
r15 = 12
go to 800
end if
************************************************************************
* *
* If there was an associated option value read from the options input *
* file, copy the option value to the local variable named optval. *
* If no associated option value was read from the options input file, *
* the option value should be set to null (i.e., a character string of *
* blanks) since for some options the option value might be blank. *
* *
************************************************************************
if (termcnt.ge.2) then
optval = optterms(2)
else
optval = ' '
end if
************************************************************************
* *
* Validate the option keyword. Once the keyword is recognized, copy *
* the associated value to the appropriate location in the option *
* common block. Then go to the top of this loop to read the next *
* line from the global options input file. *
* *
************************************************************************
if (keyword.eq.'modelname') then
modelname = optval
go to 200
end if
if (keyword.eq.'fixdname') then
fixdname = optval
go to 200
end if
if (keyword.eq.'randname') then
randname = optval
go to 200
end if
if (keyword.eq.'dumpname') then
dumpname = optval
go to 200
end if
if (keyword.eq.'logitname') then
logitname = optval
go to 200
end if
if (keyword.eq.'gibbsname') then
gibbsname = optval
go to 200
end if
if (keyword.eq.'residname') then
residname = optval
go to 200
end if
if (keyword.eq.'mlfrprioname') then
mlfrprioname = optval
go to 200
end if
if (keyword.eq.'printopt') then
printopt = optval
go to 200
end if
if (keyword.eq.'provset') then
provset = optval
go to 200
end if
************************************************************************
* *
* Continue validating the option keyword. *
* *
************************************************************************
if (keyword.eq.'fixedmethod') then
fixedmethod = optval
go to 200
end if
if (keyword.eq.'fixedprior') then
fixedprior = optval
go to 200
end if
if (keyword.eq.'fixedsupport') then
fixedsupport = optval
go to 200
end if
if (keyword.eq.'randprior') then
randprior = optval
go to 200
end if
if (keyword.eq.'varprior') then
varprior = optval
go to 200
end if
************************************************************************
* *
* The remaining options cannot accept a null (or blank) string for *
* their associated value. So if the option value is null (or blank), *
* skip the rest of the keyword validation checks and go immediately *
* to handle the error situtation. *
* *
************************************************************************
if (optval.eq.' ') go to 600
************************************************************************
* *
* If the keyword matched none of the preceding option names, then the *
* associated value needs to be checked further to see if it is valid *
* for whatever option name was specified. For example, the next set *
* of options must take on real values. To see if this might be the *
* case, try to convert the associated option value into a double *
* precision real form. If the option value can be read using a *
* double precision real format, then check to see if one of the real *
* keywords was specified. If the keyword is recognized, then the *
* double precision real option can be set according to the user's *
* input. *
* *
************************************************************************
read (optval,'(f16.0)',err=400) realval
if (keyword.eq.'fixedmargvar') then
fixedmargvar = realval
go to 200
end if
if (keyword.eq.'varshape') then
varshape = realval
go to 200
end if
if (keyword.eq.'varscale') then
varscale = realval
go to 200
end if
if (keyword.eq.'firstvariance') then
firstvariance = realval
go to 200
end if
************************************************************************
* *
* The next set of options must take on integer values. To see if *
* this might be the case, try to convert the associated option value *
* into an internal integer form. If the option value can be read *
* using an integer format, then check to see if one of the integer *
* option keywords was specified. If the keyword is recognized, then *
* the integer option can be set according to the user's input. *
* *
************************************************************************
read (optval,'(i8)',err=400) intval
if (keyword.eq.'iteracnt') then
iteracnt = intval
go to 200
end if
if (keyword.eq.'burnin') then
burnin = intval
go to 200
end if
************************************************************************
* *
* Continue validating the option keyword. *
* *
************************************************************************
if (keyword.eq.'priosize') then
priosize = intval
go to 200
end if
if (keyword.eq.'parmdps') then
parmdps = intval
go to 200
end if
if (keyword.eq.'cseed') then
cseed = intval
go to 200
end if
if (keyword.eq.'tseed') then
tseed = intval
go to 200
end if
if (keyword.eq.'randtracecnt') then
randtracecnt = intval
go to 200
end if
if (keyword.eq.'fixedcnt') then
fixedcnt = intval
go to 200
end if
if (keyword.eq.'randcnt') then
randcnt = intval
go to 200
end if
if (keyword.eq.'fixedndx') then
fixedndx = intval
go to 200
end if
if (keyword.eq.'womanid') then
womanid = intval
go to 200
end if
************************************************************************
* *
* The rest of the options are of the the 'T'rue/'F'alse variety, so *
* the first character of optval must be either a 'T' or a 'F'. *
* If this is not the case, the associated option value is invalid. *
* *
************************************************************************
400 torf = optval(1:1)
if (torf.eq.'T' .or. torf.eq.'F') then
if (keyword.eq.'ieee') then
ieee = torf
go to 200
end if
if (keyword.eq.'standardize') then
standardize = torf
go to 200
end if
if (keyword.eq.'intcdump') then
intcdump = torf
go to 200
end if
if (keyword.eq.'residout') then
residout = torf
go to 200
end if
if (keyword.eq.'gibbstrace') then
gibbstrace = torf
go to 200
end if
if (keyword.eq.'gibbstime') then
gibbstime = torf
go to 200
end if
if (keyword.eq.'mlfrpriotrace') then
mlfrpriotrace = torf
go to 200
end if
if (keyword.eq.'fixedchk') then
fixedchk = torf
go to 200
end if
if (keyword.eq.'randchk') then
randchk = torf
go to 200
end if
************************************************************************
* *
* Continue validating the option keyword. *
* *
************************************************************************
if (keyword.eq.'varchk') then
varchk = torf
go to 200
end if
if (keyword.eq.'keep4707') then
keep4707 = torf
go to 200
end if
if (keyword.eq.'longeropen') then
longeropen = torf
go to 200
end if
else
************************************************************************
* *
* If this point is reached in the program, the associated option *
* value was not valid for whatever option was named in the keyword. *
* For example, something other than a 'T' or 'F' was found where a *
* 'T'rue/'F'alse indicator was expected. In this case, set the error *
* return code to 20 before returning to the caller. *
* *
************************************************************************
600 r15 = 20
go to 800
end if
************************************************************************
* *
* If the program falls through to this point, the option keyword has *
* not been recognized. In this case, set the error return code to *
* 16 before returning to the caller. *
* *
************************************************************************
r15 = 16
go to 800
************************************************************************
* *
* If instead we have encountered an end-of-file reading the global *
* options input file set the error return code to indicate that all *
* went as would be expected before making the return to our caller. *
* *
************************************************************************
700 r15 = 0
800 close (unit=20)
return
end
************************************************************************
* *
* oneparse 2-17-94 *
* *
* This subroutine is a very rudimentary parser. The input character *
* string passed as the first argument is parsed into 0 or more tokens *
* where the tokens are separated by instances of the single character *
* delimiter passed as the delimit argument. The separate tokens are *
* returned in the output character vector, tokens. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* instring = a character string to be parsed into separate tokens. *
* *
* delimit = a single character to be used to separate the input *
* into individual tokens. *
* *
* maxtok = an integer containing the maximum number of tokens *
* which may be returned. This is the number of elements *
* available in the tokens vector. *
* *
* *
* Outputs: *
* *
* tokens = a vector of character strings containing the parsed *
* version of the input character string, instring. Each *
* entry in tokens will be left justified (leading blanks *
* removed). The number of tokens found in the input *
* string is returned in tokcnt. *
* *
* tokcnt = an integer containing the number of tokens actually *
* returned in tokens. This will be a number between 0 *
* and maxtok, inclusive. *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = the input character string contained more than *
* maxtok tokens. The first maxtok tokens have *
* been returned in tokens and tokcnt has been set *
* equal to maxtok. *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine oneparse(instring,delimit,maxtok,tokens,tokcnt,r15)
character instring*(*)
character delimit*1
integer maxtok
character tokens(maxtok)*(*)
integer tokcnt
integer r15
************************************************************************
* *
* The following variables hold various temporary values used in this *
* subroutine such as loop counters, indices and so forth. *
* *
************************************************************************
integer inlen
integer index
integer bpos
integer dpos
integer epos
integer len
integer tn
integer pn
************************************************************************
* *
* Parse the input character string into its separate tokens. The *
* input string is assumed to contain from 0 to maxtok tokens *
* separated by the token separator character, delimit. *
* *
************************************************************************
inlen = len(instring)
bpos = 1
tn = 0
************************************************************************
* *
* Locate the beginning and the end of the next token within the *
* input string. *
* *
* First, find the next nonblank character in the input string. *
* *
************************************************************************
400 do 450 pn=bpos,inlen
if (instring(pn:pn).ne.' ') go to 500
450 continue
************************************************************************
* *
* If no nonblank characters were found in the input string, there *
* are no more tokens in the input, so we have finished parsing this *
* input string. *
* *
* At this point set tokcnt to the number of tokens found and set the *
* error return code to indicate that all went as would be expected *
* before returning to the caller. *
* *
************************************************************************
tokcnt = tn
r15 = 0
return
************************************************************************
* *
* The beginning of the next token has been located. Increment the *
* subscript to use for this token. The element of the tokens vector *
* addressed by this subscript will be the one used to return the *
* current token. *
* *
* If this subscript is greater than maxtok, the passed input string *
* contained more tokens than can be returned in the tokens vector. *
* In this case, tokcnt is set equal to maxtok and the error return *
* code is set to 4 before returning to the caller. *
* *
* To find the end of this token, find the next token separator or the *
* end of the input string if there are no more separators. *
* *
************************************************************************
500 tn = tn + 1
if (tn.gt.maxtok) then
tokcnt = maxtok
r15 = 4
return
end if
bpos = pn
dpos = index(instring(bpos:),delimit)
if (dpos.eq.0) dpos = inlen
epos = bpos + dpos - 2
************************************************************************
* *
* We now have the positions of both the beginning and the end of the *
* next token within the input string, so save this token in the *
* correct element of the tokens vector. *
* *
************************************************************************
tokens(tn) = instring(bpos:epos)
************************************************************************
* *
* Start looking for the next token immediately following the most *
* recently located token separator (the one which terminated the *
* previous token). *
* *
************************************************************************
bpos = epos + 2
go to 400
end
************************************************************************
* *
* translat 11-15-92 *
* *
* Using the passed input and output translate tables, the translat *
* subroutine copies an input vector of character strings to an output *
* vector changing every instance of any character located in the *
* input translate table to the character located in the corresponding *
* position in the output table. Any character not located in the *
* input table is copied without any modification. The lengths of the *
* input and output tables must be identical. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* intable = an arbitrary length character string containing a list *
* of characters to be translated to some character. *
* Note that if the same character occurs more than once *
* in this character string, only the first instance will *
* actually be used. *
* *
* outtable = a character string of the same length as the intable *
* character string. Every time a character in the *
* intable character string is encountered in any string *
* in the vector of input character strings, instrng, it *
* is replaced (in the copy made to the output string *
* vector, outstrng) by the character in outtable which *
* is located at the position corresponding to the *
* matched character in the intable character string. *
* *
* instrng = a vector of character strings which are to be copied *
* to the output string vector, outstrng. The copy in *
* outstrng will have all instances of any character *
* located in the intable translated to the correspond- *
* ingly located character in outtable. *
* *
* strngcnt = an integer containing the actual number of character *
* strings in the instrng vector which are to be trans- *
* lated. There must be at least this many elements *
* available in the output character string vector, *
* outstrng, below. If not, this subroutine will copy *
* character strings past the end of the outstrng vector, *
* over-writing whatever is there. *
* *
************************************************************************
************************************************************************
* *
* Outputs: *
* *
* outstrng = a vector of character strings in which the input *
* vector of strings, instrng, has been copied. All *
* instances of characters found in instrng matching *
* any character in intable have been translated to the *
* correspondingly positioned character of outtable. *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = the lengths of the intable and outtable are *
* not equal. *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine translat(intable,outtable,instrng,strngcnt,outstrng,
+ r15)
character intable*(*)
character outtable*(*)
character instrng(*)*(*)
integer strngcnt
character outstrng(*)*(*)
integer r15
************************************************************************
* *
* The following variables hold various temporary values used in this *
* subroutine. *
* *
* curchar - character currently being copied (and translated) to *
* the outstrng vector *
* tablen - length of the intable and outtable character strings *
* sn - subscript of the current character string *
* cn - position of the current character (being copied) *
* tn - current position in the intable character string *
* *
************************************************************************
character curchar*1
integer tablen
integer sn
integer cn
integer tn
************************************************************************
* *
* First, check that the lengths of the intable and outtable character *
* strings are equal. *
* *
************************************************************************
tablen = len(intable)
if (len(outtable).ne.tablen) then
r15 = 4
return
end if
************************************************************************
* *
* Now perform the copying and translation of the vector of input *
* strings, instrng, into the vector of output strings, outstrng. *
* *
************************************************************************
do 600 sn=1,strngcnt
outstrng(sn) = instrng(sn)
do 500 cn=1,len(instrng(sn))
curchar = instrng(sn)(cn:cn)
do 300 tn=1,tablen
if (intable(tn:tn).eq.curchar) then
outstrng(sn)(cn:cn) = outtable(tn:tn)
go to 500
end if
300 continue
500 continue
600 continue
r15 = 0
return
end
************************************************************************
* *
* modelin 2-12-94 *
* *
* This subroutine reads the names of all covariates to be included in *
* the model from the file named as the first argument. This file is *
* referred to as the model input file. The covariate names are *
* returned in parsed representation in the character matrix passed as *
* the fourth argument and in the originally input representation in *
* the character vector passed as the fifth argument. *
* *
* The first line of the model input file is assumed to contain a left *
* justified character string of length 2 where the first character is *
* one of the integers 0 through 9 and the second character is either *
* a blank or a '+' symbol. This string indicates the caller's choice *
* of parities to include in the model. The first character integer *
* value indicates the first (chronologically earliest) parity to *
* include in the model. If the second character is a blank, then *
* only birth intervals with the parity indicated by the first *
* character are to be included in the model. If the second character *
* is a '+', then all birth intervals with parities greater than or *
* equal to the parity indicated by the first character are to be *
* included in the model. The modelin subroutine returns both parts *
* of this parity selection character string in separate output *
* arguments. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* filename = a character string containing the name of a file which *
* contains the names of the covariates to be included in *
* this run of the model. *
* *
* maxcovar = an integer containing the maximum number of covariate *
* names which can fit in the covterms matrix and in the *
* covnames vector. *
* *
* maxorder = an integer containing the highest order interaction *
* covariate permitted in the model. This is the number *
* of columns available in the covterms matrix. *
* *
************************************************************************
************************************************************************
* *
* Outputs: *
* *
* covterms = a matrix of character strings containing the parsed *
* version of the covariate variable names. Main effect *
* covariate names will be copied to the first column of *
* covterms. Interaction covariates will have the *
* individual terms separated out, one per column, up to *
* a maximum of maxorder terms per covariate. All of the *
* entries in covterms will be left justified (leading *
* blanks removed) and translated into lowercase letters. *
* The number of meaningful terms in each row (i.e., the *
* order of the covariate) is returned in the correspond- *
* ing element of the termcnt vector. *
* *
* covnames = a vector of character strings which is to be setup by *
* the modelin subroutine. This vector must contain at *
* least maxcovar+1 elements. The need to check for an *
* end-of-file explains why covnames contains an extra *
* "apparently unnecessary" appended element. The number *
* of elements actually read in is returned in the *
* numcovar argument. This vector will contain the *
* covariate identifiers in precisely the format read *
* from the model input file. No parsing or translation *
* will have been performed on this vector. If any line *
* of the model input file is longer than the covnames *
* string length, it will be truncated on the right. *
* *
* numcovar = an integer which will contain the actual number of *
* covariates the user has indicated are to be included *
* in the model. That is, numcovar will contain one less *
* than the number of lines actually read from the model *
* input file (the first line does not contain the name *
* of a covariate). *
* *
* termcnt = an integer vector containing the number of meaningful *
* terms in each row of the covterms matrix. All entries *
* in this vector will be values between 1 and maxorder, *
* inclusive. If the entry is a 1, the corresponding *
* covariate is a main effect. If it is a 2, the *
* corresponding covariate is a 2-way interaction. And *
* for a 3, the corresponding covariate is a 3-way *
* interaction, and so on. There are only numcovar valid *
* entries in this vector. *
* *
* parint = an integer containing the first (chronologically *
* earliest) parity to be included in the model. *
* *
* partype = a one character character string which will be set to *
* either a blank or a '+' character depending on which *
* is read from the model input file. *
* *
************************************************************************
************************************************************************
* *
* Outputs (continued): *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = unable to open the model input file named in *
* the filename argument. *
* 8 = an invalid parity selection argument was read *
* as the first line of the model input file. *
* 12 = more covariate names were read from the model *
* input file than could be held in the covnames *
* vector. *
* 16 = the translat subroutine returned with a nonzero *
* error return code. *
* 20 = the oneparse subroutine returned with a nonzero *
* error return code. *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine modelin(filename,maxcovar,maxorder,covterms,covnames,
+ numcovar,termcnt,parint,partype,r15)
character filename*(*)
integer maxcovar
integer maxorder
character covterms(maxcovar,maxorder)*(*)
character covnames(maxcovar+1)*(*)
integer numcovar
integer termcnt(maxcovar)
integer parint
character partype*1
integer r15
************************************************************************
* *
* The following variables hold various constants and temporary values *
* used in this subroutine such as loop counters and so forth. *
* *
* The length of the character strings in curterms should be the same *
* length as the character strings in the covterms matrix. Unfortu- *
* nately fortran does not permit dynamic memory allocation so these *
* lengths must be hard-coded at compile time. *
* *
************************************************************************
character ucase*26 /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
character lcase*26 /'abcdefghijklmnopqrstuvwxyz'/
character delimit /'*'/
character curterms(3)*18
character covcurr*80
character parin*1
integer termmax /3/
integer cn
integer tn
integer rc
************************************************************************
* *
* Open the model input file. If unable to do so, return to the *
* caller with an error code of 4. *
* *
************************************************************************
open (unit=21,file=filename,status='old',iostat=rc)
if (rc.ne.0) then
write (0,*) 'unable to open the model input file'
r15 = 4
return
end if
************************************************************************
* *
* Read the first line of the model input file as two single character *
* variables. After verifying that the second character is a valid *
* parity selection character (a blank or a '+'), convert the first *
* character into its integer form so that it may be returned to the *
* caller. *
* *
************************************************************************
read (21,'(2a1)') parin, partype
if ((partype.eq.'+' .or. partype.eq.' ') .and. parin.ge.'0' .and.
+ parin.le.'9') then
read (parin,'(i1)') parint
else
write (0,*)
+ 'invalid parity selection argument found in model file'
r15 = 8
go to 800
end if
************************************************************************
* *
* Next, read the names of the covariates to be included in this run *
* of the model. These names are stored in covnames. If there are *
* more covariates named in the model input file than the number of *
* elements in covnames, then the caller has tried to include more *
* covariates than the current implementation can handle. When this *
* occurs, modelin returns to the caller with an error return code of *
* 12. This situation is recognizable if the do-loop below ends *
* normally, that is, if no end-of-file is encountered. The need to *
* check for an end-of-file explains why the loop is executed up to *
* 'maxcovar + 1' times and why covnames contains an "apparently *
* unnecessary" element. *
* *
************************************************************************
do 200 cn=1,(maxcovar+1)
read (21,'(a)',end=300) covnames(cn)
200 continue
* As explained, it is an error if the loop has ended normally:
write (0,*)
+ 'maximum number of allowed covariates has been exceeded'
r15 = 12
go to 800
* Make sure next element of covnames won't match any covariate name.
300 covnames(cn) = ' '
numcovar = cn - 1
************************************************************************
* *
* Now parse each covariate identifier into its separate terms. The *
* delimit variable contains the special character acting as the term *
* separator. Each covariate identifier consists of from 1 to *
* maxorder terms separted by the term separator. *
* *
* While parsing the covariate identifier, also translate all terms *
* into lowercase letters. In this fashion the program will no longer *
* have case sensitive covariate identifications. *
* *
************************************************************************
do 600 cn=1,numcovar
call translat(ucase,lcase,covnames(cn),1,covcurr,rc)
if (rc.ne.0) then
write (0,*) 'translat exited with a nonzero error code of', rc
r15 = 16
go to 800
end if
call oneparse(covcurr,delimit,termmax,curterms,termcnt(cn),rc)
if (rc.ne.0) then
write (0,*) 'oneparse exited with a nonzero error code of', rc
r15 = 20
go to 800
end if
************************************************************************
* *
* Copy the resulting parsed version of the covariate identifier into *
* the vector passed to this subroutine. This copy is necessitated *
* because fortran stores matrices in column major order so that the *
* elements of each row do not occupy consecutive locations in memory. *
* *
************************************************************************
do 500 tn=1,termcnt(cn)
covterms(cn,tn) = curterms(tn)
500 continue
600 continue
************************************************************************
* *
* At this point set the error return code to indicate that all went *
* as would be expected. *
* *
************************************************************************
r15 = 0
800 close (unit=21)
return
end
************************************************************************
* *
* prvparse 4-12-93 *
* *
* This subroutine converts a character string in the format of the *
* province inclusion option located in the options common block, *
* provset, into the corresponding vector of 'T' or 'F' character *
* strings, provuse, which is used by the chkevent subroutine to *
* determine whether or not the event passed to it should be included *
* while estimating the current model. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* provset = a character string specifying which of the Iranian *
* provinces are to be included during the estimation of *
* the current model. The precise format of this string *
* is documented in the optioncb.h include file. This *
* argument to prvparse will usually be the string found *
* in the options common block, but I prefer to have it *
* passed as an argument to permit a little more *
* flexibility in using this subroutine. *
* *
* provmax = an integer containing the maximum valid Iranian *
* province number; presumably this number is 22. *
* *
* *
* Outputs: *
* *
* provuse = a vector of character strings all of length one, where *
* each element is to be set to either a 'T' (to include *
* the corresponding province) or a 'F' (to not inlcude *
* the province). The provuse character string vector *
* must consist of elements numbered starting with 0 and *
* running through at least provmax, corresponding to *
* provinces 0 through provmax. *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = the provset argument contains one or more *
* invalid characters. *
* 8 = the oneparse subroutine returned with a nonzero *
* error return code. *
* 12 = at least one province number was found to be *
* outside the valid range, 0-provmax. *
* 16 = at least one lower province limit was greater *
* than its corresponding upper province limit. *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine prvparse(provset,provmax,provuse,r15)
character provset*(*)
integer provmax
character provuse(0:provmax)*1
integer r15
************************************************************************
* *
* The following variables hold various temporary values used in this *
* subroutine such as loop counters, indices and so forth. *
* *
* rangesep - a single character to be used to separate the lower *
* and upper province limits in a single range specifier *
* delimit - a single character to be used as the separator between *
* province range specifiers *
* ranges - vector containing separated range specifiers (tokens) *
* limits - lower and upper province limit in current range token *
* *
* limmax - maximum number of province limits in one range token *
* rangemax - maximum number of range specifiers allowed in provset *
* rangecnt - number of separate range specifiers found in provset *
* limcnt - number of limits (1 or 2) in current range specifier *
* lower - lower province limit in current range specifier *
* upper - upper province limit in current range specifier *
* *
************************************************************************
character validchars*13 /'0123456789,- '/
character rangesep*1 /'-'/
character delimit*1 /','/
character ranges(23)*8
character limits(2)*2
integer limmax /2/
integer rangemax
integer rangecnt
integer limcnt
integer lower
integer upper
integer index
integer len
integer cn
integer pn
integer rn
integer rc
************************************************************************
* *
* Check that the input province selection argument does not contain *
* any invalid characters. It may only contain the ten digits, the *
* range token delimiter, the within range separator and any number of *
* blank characters. If this requirement is not satisfied, return to *
* the caller with an error return code of 4. *
* *
************************************************************************
do 100 cn=1,len(provset)
rc = index(validchars,provset(cn:cn))
if (rc.eq.0) then
write (0,*) 'provset option contains invalid characters'
r15 = 4
return
end if
100 continue
************************************************************************
* *
* Call the oneparse subroutine to separate the province selection *
* argument into range specifiers. *
* *
************************************************************************
rangemax = provmax + 1
call oneparse(provset,delimit,rangemax,ranges,rangecnt,rc)
if (rc.ne.0) then
write (0,*) 'oneparse exited with a nonzero error code of', rc
r15 = 8
return
end if
************************************************************************
* *
* Initialize the provuse argument to not include any province at all. *
* *
************************************************************************
do 200 pn=0,provmax
provuse(pn) = 'F'
200 continue
************************************************************************
* *
* Now parse each range specifier into either 1 or 2 ascii format *
* integer province limits. *
* *
************************************************************************
do 600 rn=1,rangecnt
call oneparse(ranges(rn),rangesep,limmax,limits,limcnt,rc)
if (rc.ne.0) then
write (0,*) 'oneparse exited with a nonzero error code of', rc
r15 = 8
return
end if
************************************************************************
* *
* Convert the limits into internal integer format. Note that all of *
* the characters in the limits matrix must, by now, consist of only *
* ascii digits or blanks. Make sure that all the province numbers *
* are between 0 and provmax, inclusive. If any province number is *
* outside this range, return to the caller with an error return code *
* of 12. Also, in the case where both a lower and upper limit were *
* specified in the same range specifier, make sure that the lower *
* limit is less than or equal to the upper limit. If not, return to *
* the caller with an error return code of 16. *
* *
************************************************************************
if (limits(1).eq.' ') go to 400
read (limits(1),'(i2)') lower
if (lower.lt.0) go to 400
if (limcnt.gt.1) then
read (limits(2),'(i2)') upper
if (lower.gt.upper) then
write (0,*) 'incorrectly ordered range in provset option'
r15 = 16
return
end if
else
upper = lower
end if
if (upper.gt.provmax) then
400 write (0,*) 'invalid province number in provset option'
r15 = 12
return
end if
************************************************************************
* *
* Now set each element of the provuse output argument corresponding *
* to any province number lying between the lower and upper limits, *
* inclusive, to a 'T' in order to include the province during the *
* current model estimation. *
* *
************************************************************************
do 500 pn=lower,upper
provuse(pn) = 'T'
500 continue
600 continue
************************************************************************
* *
* Everything has gone as expected if we made it to this point in the *
* program, so return to the caller with the good news. *
* *
************************************************************************
r15 = 0
return
end
************************************************************************
* *
* provin 1-23-94 *
* *
* This subroutine contains instructions to load up the province *
* character string vector, province, with a complete set of province *
* level variables for each of the 23 Iranian provinces. Each record *
* of the province.output file is read in turn. It is assumed that *
* the records in this file have already been sorted into increasing *
* province number order, from 0 to provmax (=22). The (provmax+1)th *
* element of the province character string vector is set to all NA's *
* (in case any woman is missing her province number). *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* provmax = an integer containing the maximum number province in *
* Iran (presumably this is 22). The province output *
* vector must consist of elements numbered starting with *
* 0 and running through at least (provmax+1). *
* *
* *
* Outputs: *
* *
* province = the completed province character string vector. It *
* currently must contain at least elements 0 through *
* (provmax+1), each of length at least 97. *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine provin(provmax,province,r15)
integer provmax
character province(0:provmax+1)*(*)
integer r15
************************************************************************
* *
* The following variables hold various constants and temporary values *
* used in this subroutine such as loop counters and so forth. *
* *
************************************************************************
character na(18)*2 /18*'NA'/
integer pn
integer f1
************************************************************************
* *
* Read the province variables into the province vector. *
* The input file is assumed to be sorted in ascending province number *
* order. *
* *
************************************************************************
open (unit=31,file='province.output',status='old')
do 200 pn=0,provmax
read (31,'(a)',end=210) province(pn)
200 continue
************************************************************************
* *
* We know that there are 3 NA's within the province level variables. *
* Two of these are from Boyrahmad for variables we are not using, so *
* we can just leave them in the province vector as is. However, we *
* are using the userrt variable, so the NA from Khuzestan for this *
* variable must be replaced by an imputed value. Using S, we did a *
* linear imputation using clinirt as an independent variable to find *
* a reasonable imputed value for the Khuzestan userrt variable. This *
* value is placed in the province vector here so we don't need to do *
* anything further to handle this imputation. *
* *
************************************************************************
province(6)(93:97) = '126.5'
************************************************************************
* *
* Place NA's in the appropriate columns in the last (extra) element *
* of the province vector in order to be able to produce the necessary *
* NA's to fill out the province portion for any women for which the *
* province number is not available (missing). *
* *
************************************************************************
210 write (province(provmax+1),'(5x,a,2x,a,10(3x,a),4x,a,3x,a,2(5x,a),
+ 3x,a,4x,a)') (na(f1), f1=1,18)
close (unit=31)
r15 = 0
return
end
************************************************************************
* *
* covarin 12-31-93 *
* *
* This subroutine contains instructions to load up the covariate *
* character string vector, covar, with a complete set of covariate *
* variables for each individual woman in the sample. Each record of *
* the covar.output file is read in turn. The zeroth element of the *
* covar vector is used to originally read the next record from the *
* covariate file. *
* *
* After each record is read in, repeated calls to ifsedit are made *
* to construct a standard set of additional covariates which are *
* appended to the end of each element of this covariate character *
* string vector. The additional covariates are currently: *
* *
* wedcat = woman's categorical completed education *
* hedcat = husband's categorical completed education *
* earlyage = woman's age at marriage short of 12 years *
* lateage = woman's age at marriage in excess of 12 years *
* cohort(6) = woman's cohort (a set of 6 indicator variables) *
* modern = has woman ever used modern contraception? *
* tradition = has woman ever used traditional contraception? *
* *
* After the additional covariates have been appended, the woman's id *
* is determined and the completed character string is then copied to *
* the appropriate covar element. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* womaxcnt = an integer containing the maximum number of women in *
* the sample. In other words, womaxcnt contains the *
* largest woman identification number used. It is the *
* upper limit for the subscript into the covariate *
* character string vector, covar. *
* *
* *
* Outputs: *
* *
* covar = the completed covariate character string vector. It *
* currently must contain at least elements 0 through *
* womaxcnt, each of length at least 140. *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine covarin(womaxcnt,covar,r15)
integer womaxcnt
character covar(0:womaxcnt)*(*)
integer r15
************************************************************************
* *
* The following variables hold various temporary values used in this *
* subroutine. This includes any do-loop counters and similar such *
* temporary subscripts and indices. *
* *
* inputs - character vector to pass inputs to ifsedit subroutine *
* result - character vector in which ifsedit returns its output *
* idnumb - character form of current woman identification number *
* *
* wn - main loop counter while reading from covariate file *
* id - internal integer form of woman identification number *
* r1 - subscript into the ifsedit subroutine result vector *
* s1 - first character in covar(0) to copy cohort result into *
* s2 - last character in covar(0) to copy cohort result into *
* *
************************************************************************
character inputs(6)*2
character result(6)*2
character idnumb*4
integer wn
integer id
integer r1
integer s1
integer s2
integer rc
************************************************************************
* *
* Read the separate woman covariate file into the covariate vector. *
* Each record is initially read into the 0th element so that its *
* woman id may be inspected in order to place the record in the *
* correct element of the vector. Also, a number of covariate *
* transformations and calculations are performed before placing an *
* updated version of the woman's covariates in the covariate vector. *
* *
************************************************************************
open (unit=32,file='covar.output',status='old')
do 300 wn=1,womaxcnt
read (32,'(a)',end=310) covar(0)
covar(0)(105:140) = ' '
************************************************************************
* *
* Now call ifsedit to transform the work and size covariates. *
* *
* Then, call ifsedit to calculate the 6-category woman and husband *
* education levels. *
* *
* Note: Unfortunately, I am not able to find a better way to make *
* references to the many substrings in covar(0) below while adhering *
* to the Fortran 77 standard. I would much prefer to give names to *
* each of the substrings (say, by an equivalence statement or record *
* statement), so that it would be easier to make additions and *
* deletions to the layout of covariates in the covar vector strings. *
* *
************************************************************************
call ifsedit('work',covar(0)(49:50),covar(0)(49:50),rc)
call ifsedit('size',covar(0)(52:53),covar(0)(52:53),rc)
inputs(1) = covar(0)(37:38)
inputs(2) = covar(0)(40:41)
call ifsedit('educ',inputs,covar(0)(106:107),rc)
inputs(1) = covar(0)(61:62)
inputs(2) = covar(0)(64:65)
call ifsedit('educ',inputs,covar(0)(109:110),rc)
************************************************************************
* *
* Next, call ifsedit to calculate 2 new covariates from the woman's *
* age at marriage: 1) the negative of the number of years before age *
* 12 at which she was married (for those married before this age), *
* and 2) the number of years after age 12 at which she was married *
* (for those who married after this age). *
* *
************************************************************************
call ifsedit('agem',covar(0)(9:10),result,rc)
covar(0)(112:113) = result(1)
covar(0)(115:116) = result(2)
************************************************************************
* *
* Next, call ifsedit to calculate which cohort this woman belongs to. *
* *
************************************************************************
call ifsedit('cohort',covar(0)(76:77),result,rc)
do 280 r1=1,6
s2 = r1 * 3 + 116
s1 = s2 - 1
covar(0)(s1:s2) = result(r1)
280 continue
************************************************************************
* *
* Next, call ifsedit to determine whether or not this woman has ever *
* used a modern form of contraception. *
* *
************************************************************************
inputs(1) = covar(0)(82:83)
inputs(2) = covar(0)(85:86)
inputs(3) = covar(0)(88:89)
inputs(4) = covar(0)(94:95)
call ifsedit('contraception',inputs,covar(0)(136:137),rc)
************************************************************************
* *
* Finally, call ifsedit to determine whether or not this woman has *
* ever used a traditional form of contraception. *
* *
************************************************************************
inputs(1) = covar(0)(91:92)
inputs(2) = covar(0)(97:98)
inputs(3) = covar(0)(100:101)
inputs(4) = covar(0)(103:104)
call ifsedit('contraception',inputs,covar(0)(139:140),rc)
************************************************************************
* *
* Place the just read covar string (with the updates applied by *
* ifsedit above) into the correct element of covar. *
* *
************************************************************************
idnumb = covar(0)(1:4)
read (idnumb,'(i4)') id
covar(id) = covar(0)
300 continue
310 close (unit=32)
r15 = 0
return
end
************************************************************************
* *
* ifsedit 2-06-94 *
* *
* This program is a collection of subprograms to perform a standard *
* set of edits on individual fields from the Iran Fertility Survey *
* database. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* option = a character string containing one of the following *
* edit options: *
* 'work' *
* 'size' *
* 'educ' *
* 'agem' *
* 'cohort' *
* 'contraception' *
* *
* inputs = an array of character strings containing the set of *
* input values needed to perform the selected edit. *
* *
* *
* Outputs: *
* *
* result = an array of character strings containing the *
* resulting standardized output(s) for the type of *
* edit indicated by the passed option. *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = unrecognizable option argument. *
* 8 = input data is missing (ie, is NA). *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine ifsedit(option,inputs,result,r15)
character option*(*)
character inputs(*)*(*)
character result(*)*(*)
integer r15
character workmap(9)*2 /'01', 3*'02', 3*'03', '04', '00'/
character sizemap(9)*2 /'04', '03', '02', '01', 5*'00'/
character educmap(26)*2 /5*'01', '02', 4*'00', 3*'03', 7*'00',
+ 5*'03', '04'/
integer adjage
integer argu1
integer argu2
integer i1
************************************************************************
* *
* If the input data is missing (ie, is NA), set the result to *
* missing (ie, NA) and set the return code to 8. Otherwise, convert *
* the first input character string into its internal integer form *
* and set the default error return code to 0. If need be, the *
* individual subroutines may replace this value with an error *
* condition. *
* *
************************************************************************
if ( inputs(1) .eq. 'NA' ) then
result(1) = 'NA'
r15 = 8
else
read (inputs(1),'(i2)') argu1
r15 = 0
end if
************************************************************************
* *
* Check the edit option to determine which edit to perform. For *
* the 'work' option or the 'size' option, simply transform the *
* corresponding covariates through a translate table: 'workmap' for *
* 'work' or 'sizemap' for size. *
* *
************************************************************************
if ( option .eq. 'work' ) then
if ( r15 .eq. 8 ) return
result(1) = workmap(argu1)
else if ( option .eq. 'size' ) then
if ( r15 .eq. 8 ) return
result(1) = sizemap(argu1)
************************************************************************
* *
* For the 'educ' option, the first input should be a binary valued *
* variable where 0 = no education or 1 = education. The second input *
* is the original survey form of completed level. It should be a *
* positive integer between 1 and 93. If it is between 1 and 26, use *
* the 'educmap' translate table to find the equivalent 6-category *
* value. If it is >26, it should be translated to category 5. *
* *
************************************************************************
else if ( option .eq. 'educ' ) then
if ( r15 .eq. 8 ) return
if ( inputs(2) .ne. 'NA' ) then
read (inputs(2),'(i2)') argu2
if ( argu2 .lt. 1 ) then
result(1) = 'NA'
else if ( argu2 .le. 26 ) then
result(1) = educmap(argu2)
else
result(1) = '05'
end if
else
result(1) = '00'
end if
************************************************************************
* *
* For the 'agem' option, calculate both the negative of the number *
* of years before age 12 at which the woman was married (otherwise, *
* set to 0) and the positive of the number of years after age 12 at *
* which the woman was married (otherwise, set to 0). Return both of *
* these new variables to the caller. *
* *
************************************************************************
else if ( option .eq. 'agem' ) then
if ( r15 .eq. 8 ) then
result(2) = 'NA'
return
end if
adjage = - idim(12,argu1)
write (result(1),'(i2)') adjage
adjage = idim(argu1,12)
write (result(2),'(i2)') adjage
************************************************************************
* *
* For the 'cohort' option, the first input should be an integer *
* valued variable with values between 1 and 7, inclusive. These *
* cohort levels are translated into 6 dummy variables, one for each *
* of cohorts 2 through 7. *
* *
************************************************************************
else if ( option .eq. 'cohort' ) then
if ( r15 .eq. 8 ) then
do 190 i1=2,6
result(i1) = 'NA'
190 continue
return
end if
do 200 i1=1,6
if ( i1 + 1 .eq. argu1 ) then
result(i1) = '01'
else
result(i1) = '00'
end if
200 continue
************************************************************************
* *
* For the 'contraception' option, the four input values must be *
* either 0, 1, or NA. If any of these are 1, the result should be *
* set to 1. Otherwise, if any of these are NA, the result should be *
* set to NA. Only if all four inputs are 0 should the result be set *
* to 0. *
* *
************************************************************************
else if ( option .eq. 'contraception' ) then
result(1) = '00'
do 250 i1=1,4
if ( inputs(i1) .eq. '01' ) then
result(1) = '01'
return
end if
if ( inputs(i1) .eq. 'NA' ) then
result(1) = 'NA'
end if
250 continue
************************************************************************
* *
* Otherwise, the option passed was not recognized, so return to our *
* caller with an appropriate error code. *
* *
************************************************************************
else
r15 = 4
end if
return
end
************************************************************************
* *
* chkevent 8-03-94 *
* *
* This subroutine checks that there is no missing data associated *
* with the event passed to it. It also parses the event as well as *
* obtaining all the covariates and province level data for the woman *
* to whom this event belongs. Additionally, chkevent determines *
* whether or not the event passed to it belongs to a woman residing *
* in a province which is to be included while estimating the current *
* model. It is in this subroutine that the data in the provincb, *
* covarcb and eventcb common blocks are set up. While doing this, *
* chkevent also calculates the beginning and ending years for the *
* input event, as well as the woman's identification number and year *
* of birth. These are returned to the caller. If for any reason the *
* event should not be used, a nonzero return code is passed back to *
* the caller. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* option = a character string specifying which checks within this *
* subroutine are to be performed. If this argument is *
* 'complete', all checks are to be performed. If this *
* argument is 'forbirth', only those checks needed by *
* the birthage or birthdur programs are to be performed. *
* *
* province = a character string vector, where each element contains *
* a complete set of province level data for each Iranian *
* province. Each element must be at least 97 characters *
* long. The province vector must consist of at least *
* elements 0 through provmax, for a total of provmax + 1 *
* elements in all (since there is a province 0). *
* *
* provmax = an integer containing the maximum number province in *
* Iran (presumably this is 22). The province character *
* string vector must consist of elements numbered *
* starting with 0 and running through at least provmax. *
* *
* provuse = a vector of character strings all of length one, where *
* each element is either set to 'T' or not. The provuse *
* character string vector must consist of elements *
* numbered starting with 0 and running through at least *
* provmax, corresponding to provinces 0 through provmax. *
* The event passed to chkevent should be included while *
* estimating the current model if and only if the *
* element of provuse corresponding to the province in *
* which the woman associated with this event resides is *
* set to 'T'. If not, the event will be excluded. *
* *
************************************************************************
************************************************************************
* *
* Inputs (continued): *
* *
* covar = a character string vector, where each element contains *
* a complete set of all covariate variables for each *
* individual woman in the sample. Each element must be *
* at least 140 characters long. The covar vector *
* contains womaxcnt + 1 elements where the first element *
* is addressed by subscript 0. *
* *
* womaxcnt = an integer containing the maximum number of women in *
* the sample. In other words, womaxcnt contains the *
* largest woman identification number used. It is the *
* upper limit for the subscript into the covariate *
* character string vector, covar. *
* *
* evenhist = a character string containing a complete set of all *
* variables directly associated with the event to be *
* checked out and parsed. In general, this character *
* string will be the most recently read input record *
* from the evenhist.output file. This string must *
* contain at least 39 characters. *
* *
* parint = an integer containing the caller's choice of parities *
* to include in the model. This integer indicates the *
* first (chronologically earliest) parity to include in *
* the model. *
* *
* partype = a single character which must be either a blank or a *
* "+". If this character is a blank, the caller only *
* wants events for parity parint to be included in the *
* model. If this character is a '+', then the caller *
* wants all events with parities greater than or equal *
* to parint to be included in the model. *
* *
* *
* Outputs: *
* *
* begyr = an integer containing the beginning year of the passed *
* event (in the Iranian calendar). *
* *
* endyr = an integer containing the ending year of the passed *
* event (in the Iranian calendar). *
* *
* id = an integer containing the woman identification number *
* of the woman to whom the passed event belongs. *
* *
* yob = an integer containing the year that the woman was *
* born (in the Iranian calendar). *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if nothing unusual was detected. In this *
* case, the separate subfields of evenhist have been *
* copied to eventcb and the associated variables in *
* covarcb have been set. *
* *
************************************************************************
************************************************************************
* *
* Outputs (continued): *
* *
* Otherwise, r15 can assume the following values: *
* 4 = the woman identification number passed in *
* evenhist is not in the valid range, 1-4934. *
* 8 = the passed event did not meet the parity *
* criteria specified by the caller. *
* 12 = there is at least one missing value among the *
* event specific variables required for later *
* estimation of the model. *
* 16 = there is at least one missing value among the *
* woman covariates required for estimation of the *
* model. *
* 20 = the birth interval represented by the evenhist *
* character string was found to be less than 0 *
* years long or to be a closed interval greater *
* than 24 years long or an open interval greater *
* than 38 years long. The 24 and 38 year figures *
* are meant to be (admittedly very conservative) *
* upper boundaries for lengths of closed and *
* open intervals, respectively. *
* 24 = the woman's age was found to be outside the *
* valid range, 9-51. *
* 28 = the woman's province number was found to be *
* outside the valid range, 0-22. *
* 32 = the province in which the woman was residing is *
* not to be included while estimating the current *
* model. *
* 36 = there is at least one missing value among the *
* province level data which is associated with *
* this woman and is also required for estimation *
* of the model. *
* *
* In all cases when r15 is not 0, the caller cannot *
* safely use any of the other output variables. One or *
* more of them may not have been set. *
* *
************************************************************************
subroutine chkevent(option,province,provmax,provuse,covar,
+ womaxcnt,evenhist,parint,partype,begyr,endyr,id,yob,r15)
character option*8
integer provmax
character province(0:provmax)*(*)
character provuse(0:provmax)*1
integer womaxcnt
character covar(0:womaxcnt)*(*)
character evenhist*(*)
integer parint
character partype*1
integer begyr
integer endyr
integer id
integer yob
integer r15
************************************************************************
* *
* The province character string vector, province, contains all of the *
* province level variables for each of the 23 Iranian provinces. It *
* is expedient to use a read statement to substring the vector *
* element containing the province level variables associated with the *
* current birth interval into its separate subfields. The outputs *
* from the read statement are stored in the character strings located *
* in the provincb common block. The data in this common block are *
* also used by the rowgen subroutine. Note that these declarations *
* are listed in the order in which the subfields exist in the *
* province character string. *
* *
************************************************************************
common /provincb/ educ1, illratio, malelit, married, mlnonag,
+ fmtoml, femsal, electr, radio, tv, phone, urban3, health1,
+ healtha, uhincome, rhincome, clinirt, userrt
character educ1*4
character illratio*3
character malelit*4
character married*4
character mlnonag*4
character fmtoml*4
character femsal*4
character electr*4
character radio*4
character tv*4
character phone*4
character urban3*4
character health1*5
character healtha*4
character uhincome*6
character rhincome*6
character clinirt*4
character userrt*5
************************************************************************
* *
* The covariate character string vector, covar, contains a complete *
* set of covariate variables for each individual woman in the *
* sample. It is expedient to use a read statement to substring the *
* vector element containing the covariates associated with the *
* current birth interval into its separate subfields. The outputs *
* from the read statement are stored in the character strings located *
* in the covarcb common block. The data in this common block are *
* also used by the rowgen subroutine. Note that these declarations *
* are listed in the order in which the subfields exist in the covar *
* character string. *
* *
************************************************************************
common /covarcb/ idnumb, yearb, agem, provnumb, hozeh, q401, q504,
+ q321, q322, q324, q325, q106, q107, q108, q109, work, size,
+ q103, q104, q709, q710, q711, q712, q714, agex, iage, pilluse,
+ iud, otherm, douche, condom, rhythm, withdraw, abstain, wedcat,
+ hedcat, earlyage, lateage, cohort, modern, tradition
character idnumb*4
character yearb*2
character agem*2
character provnumb*2
character hozeh*3
character q401*2
character q504*2
character q321*2
character q322*2
character q324*2
character q325*2
character q106*2
character q107*2
character q108*2
character q109*2
character work*2
character size*2
character q103*2
character q104*2
character q709*2
character q710*2
character q711*2
character q712*2
character q714*2
character agex*2
character iage*2
character pilluse*2
character iud*2
character otherm*2
character douche*2
character condom*2
character rhythm*2
character withdraw*2
character abstain*2
character wedcat*2
character hedcat*2
character earlyage*2
character lateage*2
character cohort(6)*2
character modern*2
character tradition*2
************************************************************************
* *
* Only a single record of the birth event history file will ever be *
* needed at any given point in this program. Therefore, no event *
* history vector need be declared. Instead, simply a scalar *
* character string, evenhist, is declared. It will be expedient to *
* use a read statement to substring the current birth event history *
* record into its separate subfields. The outputs from the read *
* statement are stored in the character strings located in the *
* eventcb common block. The data in this common block are also used *
* by the rowgen subroutine. Note that these declarations are listed *
* in the order in which the subfields exist in the evenhist *
* character string. *
* *
************************************************************************
common /eventcb/ intbegcm, intbegyr, intendcm, intendyr, closed,
+ parity, previous, prevlivm, prevlivy, breastfd, brestfln
character intbegcm*3
character intbegyr*2
character intendcm*3
character intendyr*2
character closed*1
character parity*2
character previous*2
character prevlivm*3
character prevlivy*2
character breastfd*2
character brestfln*2
************************************************************************
* *
* The fortran commands defining the option common block are in a *
* separate include file, 'optioncb.h'. Use of this include file is *
* not in accordance with the Fortran 77 Ansi standard. *
* *
************************************************************************
include 'optioncb.h'
************************************************************************
* *
* The following variables hold various temporary values used in this *
* subroutine. This includes any do-loop counters and similar such *
* temporary subscripts and indices. *
* *
************************************************************************
integer lenyr
integer pari
integer c1
integer pn
************************************************************************
* *
* First parse the event passed as evenhist into all of its subfields. *
* Skip this event if any of the required fields are missing. *
* *
************************************************************************
read (evenhist,'(a,11(1x,a))') idnumb, intbegcm, intbegyr,
+ intendcm, intendyr, closed, parity, previous, prevlivm,
+ prevlivy, breastfd, brestfln
if (parity.eq.'NA' .or. intbegyr.eq.'NA' .or.
+ intendyr.eq.'NA') then
r15 = 12
return
end if
************************************************************************
* *
* We only have to validate the parity and previous variables when the *
* option argument is 'complete'. *
* *
* The variable previous has no meaning for parity 0 intervals. *
* *
************************************************************************
read (parity,'(i2)') pari
if (option.eq.'complete' .and. pari.ne.0 .and. previous.eq.'NA')
+ then
r15 = 12
return
end if
************************************************************************
* *
* Determine whether or not this birth interval satisfies the user's *
* choice of parities to be included. If not, return to the caller *
* with an error code of 8. *
* *
************************************************************************
if ( (partype.eq.'+' .and. pari.lt.parint) .or.
+ (partype.eq.' ' .and. pari.ne.parint) ) then
r15 = 8
return
end if
************************************************************************
* *
* Convert character formatted intbegyr and intendyr into internal *
* integer format and calculate length of interval (in years). Do *
* not include this birth interval if it is less than 0 years long *
* or if it is a closed interval greater than 24 years long or if it *
* is an open interval greater than 38 years long. The 24 and 38 *
* year figures are meant to be rather conservative conceivable upper *
* boundaries. *
* For the time being I am not distinguishing between open and closed *
* intervals (for the Demand vs. Ideation paper). Eventually this *
* should be corrected. *
* *
************************************************************************
read (intbegyr,'(i2)') begyr
read (intendyr,'(i2)') endyr
lenyr = endyr - begyr
if (longeropen.eq.'T') then
if (lenyr.lt.0 .or. (closed.eq.'1' .and. lenyr.gt.24) .or.
+ (closed.eq.'0' .and. lenyr.gt.38)) then
r15 = 20
return
end if
else
if (lenyr.lt.0 .or. lenyr.gt.24) then
r15 = 20
return
end if
end if
************************************************************************
* *
* Obtain the covariates to be used in the model for the woman whose *
* birth interval is being processed. We might as well parse this *
* woman's covar character string into all of its separate covariates *
* and then use only those covariates we wish to keep in the model. *
* *
* If the woman id found in evenhist is outside the valid range, then *
* return to our caller with error return code 4. *
* *
* Since the logistic regression subroutine does not accept missing *
* data (NA's) in the model covariates, we need to bypass any interval *
* missing one or more of the possibly required covariates. *
* *
************************************************************************
if ('0001'.le.idnumb .and. idnumb.le.'4934') then
read (idnumb,'(i4)') id
else
r15 = 4
return
end if
read (covar(id),'(4x,45(1x,a))') yearb, agem, provnumb, hozeh,
+ q401, q504, q321, q322, q324, q325, q106, q107, q108, q109,
+ work, size, q103, q104, q709, q710, q711, q712, q714, agex,
+ iage, pilluse, iud, otherm, douche, condom, rhythm, withdraw,
+ abstain, wedcat, hedcat, earlyage, lateage, (cohort(c1),c1=1,6),
+ modern, tradition
************************************************************************
* *
* The birthage and birthdur subroutines do not refer to as many of *
* the woman associated covariates, so it is not necessary to check *
* for missing values of those covariates not required. Otherwise, *
* when the option argument is 'complete', it is necessary to check *
* to see if any of these additional covariates is missing. *
* *
* Actually birthage and birthdur use slightly different sets of *
* covariates, so the following still requires a little more work. *
* *
************************************************************************
if (yearb.eq.'NA' .or. agem.eq.'NA' .or. agex.eq.'NA') then
r15 = 16
return
end if
if (option.eq.'complete' .and. (size.eq.'NA' .or. q104.eq.'NA'
+ .or. wedcat.eq.'NA' .or. hedcat.eq.'NA' .or. modern.eq.'NA'))
+ then
r15 = 16
return
end if
************************************************************************
* *
* Convert character formatted yearb into internal integer format. *
* Don't include this birth interval if either the beginning or ending *
* age is negative (we already know that ending age is greater than or *
* equal to beginning age), or if ending age is greater than 56. *
* *
************************************************************************
read (yearb,'(i2)') yob
if (begyr.lt.yob .or. (endyr-yob).gt.56) then
r15 = 24
return
end if
************************************************************************
* *
* Obtain the province level data which is associated with this woman *
* and is also required for estimation of the model. We might as well *
* parse the associated province character string into all of its *
* separate data fields and then use only those province variables we *
* wish to incorporate in the model. *
* *
* If the province number found in this woman's covariate character *
* string is not a valid subscript into the province character string *
* vector, then return to our caller with error return code 28. *
* *
************************************************************************
if (option.eq.'complete') then
if ('00'.le.provnumb .and. provnumb.le.'22') then
read (provnumb,'(i2)') pn
else
if (keep4707.eq.'T' .and. id.eq.4707) then
pn = 16
else
r15 = 28
return
end if
end if
************************************************************************
* *
* If this woman resides in a province which is not to be included for *
* estimation of the current model, then return to our caller with *
* error return code 32. *
* *
************************************************************************
if (provuse(pn).ne.'T') then
r15 =32
return
end if
read (province(pn),'(2x,18(1x,a))') educ1, illratio, malelit,
+ married, mlnonag, fmtoml, femsal, electr, radio, tv, phone,
+ urban3, health1, healtha, uhincome, rhincome, clinirt, userrt
************************************************************************
* *
* Since the logistic regression and Metropolis procedure subroutines *
* do not accept missing data (NA's) in the model province level *
* variables, we need to bypass any interval missing one or more of *
* the potentially required province level variables. *
* *
************************************************************************
if (educ1.eq.'NA') then
r15 = 36
return
end if
end if
************************************************************************
* *
* At this point in the subroutine all the requisite data associated *
* with the event record passed to it has been validated and copied to *
* the appropriate common block. So return to our caller with return *
* code 0 (indicating that no errors were found). *
* *
************************************************************************
r15 = 0
return
end
************************************************************************
* *
* rowgen 2-07-94 *
* *
* The sole function of this subroutine is to correctly calculate and *
* store the entire set of covariates for the specified row of the *
* large internal event history covariate array, intcovar. It uses *
* inputs from three common blocks, provincb, covarcb and eventcb, *
* along with its arguments to perform this task. Note that this is a *
* very trusting subroutine. It does almost no error checking. It is *
* the responsibility of the caller to make sure that subscripts into *
* intcovar are not out of bounds. It is also the caller's respon- *
* sibility to ensure that none of the input data values required to *
* calculate any of the covariates named in the covterms matrix are *
* missing, i.e., that none of them are NA. This subroutine does not *
* check for NA's. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* covterms = a matrix of character strings containing the parsed *
* version of the covariate variable names. There are *
* maxcovar rows in covterms where each row represents *
* a covariate to be included in this run of the model. *
* Main effect covariate identifiers are passed in the *
* first column of covterms. Interaction covariates have *
* the individual terms separated out, one per column. *
* All entries in covterms should be left justified *
* (leading blanks removed) and in lowercase letters. *
* The number of meaningful rows in covterms is passed in *
* numcovar. The number of meaningful terms in each row *
* (i.e., the order of the covariate) is passed in the *
* corresponding element of the termcnt vector. *
* *
* maxcovar = an integer containing the maximum number of covariate *
* identifiers which can fit in the covterms matrix. *
* This is total number of rows available in covterms. *
* *
* termcnt = an integer vector containing the number of meaningful *
* terms in each row of the covterms matrix. All entries *
* in this vector must be values between 1 and the number *
* of columns in covterms, inclusive. If the entry is a *
* 1, the corresponding covariate is a main effect. If *
* it is a 2, the corresponding covariate is a 2-way *
* interaction. And for a 3, the corresponding covariate *
* is a 3-way interaction, and so on. There are only *
* numcovar meaningful entries in this vector. *
* *
* numcovar = an integer containing the number of covariates to *
* include in this run of the model. In other words, *
* numcovar is the number of meaningful rows in covterms *
* as well as the number of columns of intcovar which *
* need to be calculated and stored by rowgen. *
* *
* maxint = an integer containing the allocated (i.e., maximum *
* available) number of rows in the covariates matrix, *
* intcovar. *
* *
************************************************************************
************************************************************************
* *
* Inputs (continued): *
* *
* icr = an integer containing the number of the row (origin 1) *
* in intcovar which rowgen is to fill in. *
* *
* duration = an integer containing the number of years since the *
* woman's most recent birth (or marriage). *
* *
* iyear = an integer containing the Iranian calendar year *
* associated with this exposure-year. The fact that the *
* Iranian year starts on March 21st is ignored. The *
* relationship between Iranian and Western calendar *
* years is discussed in more detail later in this *
* subroutine. *
* *
* age = an integer containing the woman's age during the *
* current exposure-year. *
* *
* *
* Outputs: *
* *
* intcovar = a double precision matrix containing the complete set *
* of covariate data for all women exposure-years in the *
* sample. This matrix should be maxint rows by numcovar *
* columns in dimension. The subscript of the row to be *
* filled in is passed as icr. *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. In this *
* case, the specified row of intcovar will have been *
* completely filled in with the covariate values for *
* the current woman exposure-year. *
* *
* Otherwise, r15 can assume the following values: *
* 4 = an unrecognizable covariate identifier was *
* found in the covterms matrix. In this case, *
* intcovar will not have been completely filled *
* in. The caller should no longer trust intcovar *
* to have meaningful data in it. *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine rowgen(covterms,maxcovar,termcnt,numcovar,maxint,icr,
+ duration,iyear,age,intcovar,r15)
character covterms(maxcovar,*)*(*)
integer maxcovar
integer termcnt(numcovar)
integer numcovar
integer maxint
integer icr
integer duration
integer iyear
integer age
double precision intcovar(maxint,numcovar)
integer r15
************************************************************************
* *
* The province character string vector, province, contains all of the *
* province level variables for each of the 23 Iranian provinces. It *
* is expedient to use a read statement to substring the vector *
* element containing the province level variables associated with the *
* current birth interval into its separate subfields. The outputs *
* from the read statement are stored in the character strings located *
* in the provincb common block. The data in this common block are *
* also used by the rowgen subroutine. Note that these declarations *
* are listed in the order in which the subfields exist in the *
* province character string. *
* *
************************************************************************
common /provincb/ educ1, illratio, malelit, married, mlnonag,
+ fmtoml, femsal, electr, radio, tv, phone, urban3, health1,
+ healtha, uhincome, rhincome, clinirt, userrt
character educ1*4
character illratio*3
character malelit*4
character married*4
character mlnonag*4
character fmtoml*4
character femsal*4
character electr*4
character radio*4
character tv*4
character phone*4
character urban3*4
character health1*5
character healtha*4
character uhincome*6
character rhincome*6
character clinirt*4
character userrt*5
************************************************************************
* *
* The covariate character string vector, covar, contains a complete *
* set of covariate variables for each individual woman in the *
* sample. It is expedient to use a read statement to substring the *
* vector element containing the covariates associated with the *
* current birth interval into its separate subfields. The outputs *
* from the read statement are stored in the character strings in the *
* covarcb common block. Note that these declarations are listed in *
* the order in which the subfields exist in the covar character *
* string. *
* *
************************************************************************
common /covarcb/ idnumb, yearb, agem, provnumb, hozeh, q401, q504,
+ q321, q322, q324, q325, q106, q107, q108, q109, work, size,
+ q103, q104, q709, q710, q711, q712, q714, agex, iage, pilluse,
+ iud, otherm, douche, condom, rhythm, withdraw, abstain, wedcat,
+ hedcat, earlyage, lateage, cohort, modern, tradition
character idnumb*4
character yearb*2
character agem*2
character provnumb*2
character hozeh*3
character q401*2
character q504*2
character q321*2
character q322*2
character q324*2
character q325*2
character q106*2
character q107*2
character q108*2
character q109*2
character work*2
character size*2
character q103*2
character q104*2
character q709*2
character q710*2
character q711*2
character q712*2
character q714*2
character agex*2
character iage*2
character pilluse*2
character iud*2
character otherm*2
character douche*2
character condom*2
character rhythm*2
character withdraw*2
character abstain*2
character wedcat*2
character hedcat*2
character earlyage*2
character lateage*2
character cohort(6)*2
character modern*2
character tradition*2
************************************************************************
* *
* Only a single record of the birth event history file will ever be *
* needed at any given point in this program. Therefore, no event *
* history vector need be declared. Instead, simply a scalar *
* character string, evenhist, is declared. It will be expedient to *
* use a read statement to substring the current birth event history *
* record into its separate subfields. The outputs from the read *
* statement are stored in the character strings declared following *
* the evenhist character string. Note that these declarations are *
* listed in the order in which the subfields exist in the evenhist *
* character string. *
* *
************************************************************************
common /eventcb/ intbegcm, intbegyr, intendcm, intendyr, closed,
+ parity, previous, prevlivm, prevlivy, breastfd, brestfln
character intbegcm*3
character intbegyr*2
character intendcm*3
character intendyr*2
character closed*1
character parity*2
character previous*2
character prevlivm*3
character prevlivy*2
character breastfd*2
character brestfln*2
************************************************************************
* *
* The compare function is a double precision function which returns *
* a 1.0d0 if a given relationship between integer arguments pertains, *
* and otherwise returns a 0.0d0. *
* *
************************************************************************
double precision compare
************************************************************************
* *
* Using the IFS dataset, we have derived a functional approximation *
* to the expected fertility rates for durations 3 and greater. The *
* dur3plus function returns this functional approximation given the *
* current value of the duration variable. The calculation of this *
* approximation is required at several different locations during the *
* calculations of the covariates in the rowgen subroutine. *
* *
************************************************************************
double precision dur3plus
************************************************************************
* *
* The usescale function multiplies whichever contraception-ever-used *
* variable is passed to it by a factor representing the relative *
* proportion of modern contraception users in the current exposure *
* year. *
* *
************************************************************************
double precision usescale
************************************************************************
* *
* The primeduc vector contains the proportion of the eligible popula- *
* tion attending primary education during each valid exposure-year. *
* *
************************************************************************
double precision primeduc(4:56) / 0.029212d0, 0.034715d0,
+ 0.040306d0, 0.045226d0, 0.050000d0, 0.054634d0, 0.060577d0,
+ 0.062085d0, 0.065421d0, 0.069585d0, 0.076364d0, 0.078853d0,
+ 0.092593d0, 0.105469d0, 0.111966d0, 0.111532d0, 0.118988d0,
+ 0.103758d0, 0.098229d0, 0.094444d0, 0.100939d0, 0.113251d0,
+ 0.130930d0, 0.144353d0, 0.199042d0, 0.275073d0, 0.232891d0,
+ 0.258304d0, 0.260021d0, 0.264443d0, 0.277231d0, 0.304552d0,
+ 0.326087d0, 0.347307d0, 0.384683d0, 0.403608d0, 0.421638d0,
+ 0.449060d0, 0.463962d0, 0.494401d0, 0.513654d0, 0.547560d0,
+ 0.564693d0, 0.583881d0, 0.598768d0, 0.597612d0, 0.623986d0,
+ 0.650926d0, 0.673937d0, 0.745386d0, 0.791918d0, 0.843819d0,
+ 0.899123d0/
************************************************************************
* *
* The femalebyprov vector contains the proportion of eligible female *
* children attending primary education in each of the 23 Iranian *
* provinces as determined during the 1976 census (provided by Akbar *
* Aghajanian in email dated May 24, 1992). *
* *
************************************************************************
double precision femalebyprov(0:22) / 0.805d0, 0.779d0, 0.685d0,
+ 0.421d0, 0.523d0, 0.577d0, 0.645d0, 0.675d0, 0.643d0, 0.516d0,
+ 0.694d0, 0.453d0, 0.318d0, 0.384d0, 0.457d0, 0.447d0, 0.360d0,
+ 0.401d0, 0.524d0, 0.334d0, 0.826d0, 0.767d0, 0.451d0/
************************************************************************
* *
* The malebyprov vector contains the proportion of eligible male *
* children attending primary education in each of the 23 Iranian *
* provinces as determined during the 1976 census (provided by Akbar *
* Aghajanian in email dated May 24, 1992). *
* *
************************************************************************
double precision malebyprov(0:22) / 0.897d0, 0.915d0, 0.884d0,
+ 0.749d0, 0.792d0, 0.836d0, 0.801d0, 0.851d0, 0.795d0, 0.775d0,
+ 0.896d0, 0.655d0, 0.762d0, 0.797d0, 0.781d0, 0.698d0, 0.653d0,
+ 0.676d0, 0.819d0, 0.715d0, 0.923d0, 0.928d0, 0.705d0/
************************************************************************
* *
* The dureigen matrix contains the eigenvectors of the covariance *
* matrix associated with the four duration covariates. The eigen- *
* vectors are stored in column format. *
* *
************************************************************************
double precision dureigen(0:3,4) / 0.72989d0, -0.68287d0,
+ -0.03003d0, -0.00651d0, -0.43017d0, -0.49353d0, 0.75303d0,
+ 0.06575d0, 0.41468d0, 0.42388d0, 0.56480d0, -0.57391d0,
+ 0.33204d0, 0.33233d0, 0.33622d0, 0.81625d0/
************************************************************************
* *
* The ethnicity vector contains integers between 0 and 5, inclusive, *
* indicating which is the majority ethnic group in each of the 23 *
* Iranian provinces as determined during the 1976 census (provided by *
* Akbar Aghajanian in email dated May 24, 1992). The mapping from *
* the integers to the ethnic groups is as follows: 0 -> arab, *
* 1 -> baluch, 2 -> kurd, 3 -> lor, 4 -> persian, 5 -> turk. *
* *
************************************************************************
integer ethnicity(0:22) /4, 4, 4, 5, 5, 2, 0, 4, 4, 4, 4, 1, 2,
+ 4, 3, 3, 2, 5, 0, 5, 4, 4, 0/
************************************************************************
* *
* The following variables hold various temporary values used in *
* different segments of this program. Centage is woman's current *
* age minus an "average age", which we chose to be 31 years old. *
* Rate is an intermediate variable for a number of calculated rates. *
* Curname contains the identifier of the individual covariate *
* currently being inserted into the design matrix of all woman-year *
* covariates, intcovar. Curlower and curupper are character repre- *
* sentations of the lower and upper limits of a range of acceptable *
* numeric values. Onechar is a single character string reserved for *
* very local temporary use when needed. These last three variables *
* are used in several of the covariate routines which follow. *
* *
************************************************************************
double precision centage
double precision curterm
double precision product
double precision rate
character curname*18
character curlower*2
character curupper*2
character twochar*2
character onechar*1
************************************************************************
* *
* Note: The years input from the IFS database are expressed in the *
* Iranian calendar. These differ from Western calendar years by 21 *
* in the last two digits of the year. For example, the Iranian year *
* 1347 roughly coincides with the Western year 1968 (actually, the *
* Iranian year starts on March 21, but this fact is ignored here). *
* *
************************************************************************
integer yearsliv
integer chngpnt
integer wyear
integer digit
integer lower
integer upper
integer isize
integer pari
integer pcn
integer pvn
integer cn
integer tn
************************************************************************
* *
* A couple variables are used by several of the following covariate *
* calculation routines. These are initialized before starting the *
* main loop below. *
* *
************************************************************************
centage = dble(age - 31) / 10.0d0
wyear = iyear + 21
************************************************************************
* *
* The set of covariates to include in the calculation of the model is *
* specified by the names in the covterms matrix. In the following *
* nested loops, each covariate identifier in covterms is examined. *
* Then, depending on the actual identifier, one (and only one) of the *
* numerous covariate calculation if-clauses is executed. I have *
* intentionally made the inclusion of each covariate its own separate *
* if-clause so that I may easily add or delete covariates from the *
* model. *
* *
************************************************************************
do 500 cn=1,numcovar
product = 1.0d0
do 400 tn=1,termcnt(cn)
curname = covterms(cn,tn)
************************************************************************
* *
* If curname is the 'intercept' identifier, include an intercept. *
* *
************************************************************************
if (curname.eq.'intercept') then
curterm = 1.0d0
go to 300
end if
************************************************************************
* *
* If curname is the 'age correction' identifier, calculate the *
* expected log fertility for a woman at this age (the coefficients *
* used here were found by performing a weighted least squares *
* regression fit of the natural logarithm of expected fertility *
* rate against centered age and centered age squared, whose values *
* were empirically estimated from the IFS dataset. Originally a *
* cubic term was also included; however it added very little to the *
* fit and so was dropped from the regression). *
* *
************************************************************************
if (curname.eq.'age correction') then
curterm = - (0.467796d0 * centage**2 +
+ 0.509182d0 * centage + 1.155003d0)
go to 300
end if
************************************************************************
* *
* If curname is the 'cubic correction' identifier, calculate the *
* expected fertility for a woman at this age (the coefficients used *
* here were found by using Cricket-Graph on a Macintosh to find a *
* best cubic polynomial fit to the sample average fertility rate *
* at this age). This covariate should be included in the model on *
* the logit scale. *
* *
************************************************************************
if (curname.eq.'cubic correction') then
rate = 2.8563d-5 * age ** 3 - 3.3419d-3 * age ** 2 +
+ 0.11202d0 * age - 0.80948d0
curterm = rate / (1.0d0 - rate)
go to 300
end if
************************************************************************
* *
* Alternatively, it might be preferable to include any or all of the *
* following three covariate functions of current age: 'cent.age', *
* 'cent.age.square', or 'adjus.age.cubed'. The middle woman's age *
* is 31. The coefficient used in calculating 'adjus.age.cubed' was *
* selected to cause these three covariates to be mutually orthogonal. *
* *
************************************************************************
if (curname.eq.'cent.age') then
curterm = centage
go to 300
end if
if (curname.eq.'cent.age.square') then
curterm = centage**2
go to 300
end if
if (curname.eq.'adjus.age.cubed') then
curterm = centage**3 - 2.518d0 * centage
go to 300
end if
************************************************************************
* *
* Alternatively, woman's age during the current exposure year may *
* be included in the model as a set of seven indicator variables. *
* In this implementation, there are separate indicators for ages *
* 15-19, ages 20-24, ages 25-29, ages 30-34, ages 35-39, ages 40-44, *
* and a last indicator for ages 45 and above. There is no indicator *
* for ages less than 15 years old. This last set of exposures is *
* being used as a baseline and as such is incorporated in the *
* intercept. *
* *
* If curname is any of the following, include the corresponding *
* indicator in the calculation of the model: *
* ages 15-19, ages 20-24, ages 25-29, ages 30-34, ages 35-39, *
* ages 40-44, ages 45+ *
* *
************************************************************************
if (curname.eq.'ages 15-19') then
curterm = compare('ge',age,15) * compare('le',age,19)
go to 300
end if
if (curname.eq.'ages 20-24') then
curterm = compare('ge',age,20) * compare('le',age,24)
go to 300
end if
if (curname.eq.'ages 25-29') then
curterm = compare('ge',age,25) * compare('le',age,29)
go to 300
end if
if (curname.eq.'ages 30-34') then
curterm = compare('ge',age,30) * compare('le',age,34)
go to 300
end if
if (curname.eq.'ages 35-39') then
curterm = compare('ge',age,35) * compare('le',age,39)
go to 300
end if
if (curname.eq.'ages 40-44') then
curterm = compare('ge',age,40) * compare('le',age,44)
go to 300
end if
if (curname.eq.'ages 45+') then
curterm = compare('ge',age,45)
go to 300
end if
************************************************************************
* *
* The expected fertility for duration of the current interval is *
* implemented using four covariates. The first three are indicator *
* variables for duration 0, duration 1, and duration 2. Durations 3 *
* and greater are treated separately in the other duration covariate. *
* We have found a functional approximation to the expected fertility *
* for durations 3 and greater. The dur3plus function returns this *
* functional approximation when passed the current value of duration. *
* *
* The four duration covariates are identified by the names: *
* duration0, duration1, duration2, and duration3+ *
* *
************************************************************************
if (curname(1:8).eq.'duration'.and.curname(11:18).eq.' ') then
if (curname(10:10).eq.' ') then
onechar = curname(9:9)
if (onechar.ge.'0' .and. onechar.le.'2') then
read (onechar,'(i1)') digit
curterm = compare('eq',duration,digit)
go to 300
end if
else
if (curname(9:10).eq.'3+') then
if (duration.ge.3) then
curterm = dur3plus(duration)
else
curterm = 0.0d0
end if
go to 300
end if
end if
end if
************************************************************************
* *
* The preceding four duration covariates turn out to be fairly highly *
* correlated with each other which leads to highly correlated *
* parameter estimates. As an alternative to including the four *
* duration covariates, models might be better fit by including one or *
* more of the principal components determined from the duration *
* covariates in the model rather than using the original variables. *
* *
* The four duration covariates principal component variables are *
* identified by the names: *
* duration pc1, duration pc2, duration pc3, and duration pc4 *
* *
************************************************************************
if (curname(1:11).eq.'duration pc' .and. curname(13:18).eq.
+ ' ') then
onechar = curname(12:12)
if (onechar.ge.'1' .and. onechar.le.'4') then
read (onechar,'(i1)') pcn
if (duration.ge.3) then
curterm = dur3plus(duration) * dureigen(3,pcn)
else
curterm = dureigen(duration,pcn)
end if
go to 300
end if
end if
************************************************************************
* *
* Parity may be included in the model using up to eleven indicator *
* variables. There are separate indicators available for each of *
* parities 0, 1, 2, 3, 4, 5, 6, 7, 8, and 9. Parities 10 and greater *
* are lumped together in one indicator variable. In each model, one *
* selected parity (usually the smallest number included in the model) *
* should be used as the baseline and as such must be incorporated in *
* the intercept. *
* *
* The eleven parity covariates available are identified by the names: *
* parity0, parity1, parity2, parity3, parity4, parity5, parity6, *
* parity7, parity8, parity9, and parity10+ *
* *
************************************************************************
if (curname(1:6).eq.'parity' .and. curname(10:18).eq.' ')
+ then
read (parity,'(i2)') pari
if (curname(7:9).eq.' ') then
curterm = dble(pari)
go to 300
end if
onechar = curname(7:7)
if (onechar.ge.'0' .and. onechar.le.'9') then
if (curname(8:9).eq.' ') then
read (onechar,'(i1)') digit
curterm = compare('eq',pari,digit)
go to 300
end if
if (curname(7:9).eq.'10+') then
curterm = compare('ge',pari,10)
go to 300
end if
end if
end if
************************************************************************
* *
* Cohort may be included in the model using six indicator variables. *
* These six indicators are identified by the names: *
* cohort 20-24, cohort 25-29, cohort 30-34, cohort 35-39, *
* cohort 40-44, and cohort 45-49. *
* The cohort of women less than 20 years old at the time of the *
* survey is being used as the baseline and as such is incorporated *
* in the intercept. *
* *
************************************************************************
if (curname.eq.'cohort 20-24') then
read (cohort(1),'(f2.0)') curterm
go to 300
end if
if (curname.eq.'cohort 25-29') then
read (cohort(2),'(f2.0)') curterm
go to 300
end if
if (curname.eq.'cohort 30-34') then
read (cohort(3),'(f2.0)') curterm
go to 300
end if
if (curname.eq.'cohort 35-39') then
read (cohort(4),'(f2.0)') curterm
go to 300
end if
if (curname.eq.'cohort 40-44') then
read (cohort(5),'(f2.0)') curterm
go to 300
end if
if (curname.eq.'cohort 45-49') then
read (cohort(6),'(f2.0)') curterm
go to 300
end if
************************************************************************
* *
* Alternatively, cohort may be included in the model as a single *
* covariate taking values 1, 2, 3, 4, 5, 6, or 7 (ie., the agex *
* covariate). This form of cohort is to be included if curname *
* is the 'cohort' identifier. *
* *
************************************************************************
if (curname.eq.'cohort') then
read (agex,'(f2.0)') curterm
go to 300
end if
************************************************************************
* *
* Period may be incorporated in the model using any collection of *
* indicator variables, subject to the constraints described below. *
* The indicator variables are identified by names of the form: *
* '19xx-19yy', where xx and yy are numbers with xx less than or *
* equal to yy. The intervals indicated by the set of xx and yy's *
* should not be overlapping. Also, it is necessary that at least *
* one period not be specified in the set of covariates so that the *
* non-included period may serve as a baseline and as such be *
* incorporated in the intercept. *
* *
************************************************************************
if (curname(1:2).eq.'19' .and. curname(5:7).eq.'-19' .and.
+ curname(10:18).eq.' ') then
curlower = curname(3:4)
curupper = curname(8:9)
if (curlower.ge.'00' .and. curlower.le.'77' .and.
+ curupper.ge.'00' .and. curupper.le.'77' .and.
+ curlower.le.curupper) then
read (curlower,'(i2)',err=220) lower
read (curupper,'(i2)',err=220) upper
curterm = compare('ge',wyear,lower) *
+ compare('le',wyear,upper)
go to 300
end if
end if
220 continue
************************************************************************
* *
* If curname is any of the following: 'pre1957 period', 'post1957 *
* period', 'pre1959 period', 'post1959 period', 'pre1961 period' or *
* 'post1961 period', include the corresponding functional form for *
* the period effect in the calculation of the model. *
* *
************************************************************************
if (curname.eq.'pre1957 period' .or. curname.eq.
+ 'pre1959 period' .or. curname.eq.'pre1961 period') then
read (curname(6:7),'(i2)') chngpnt
if (wyear.le.chngpnt) then
curterm = dble(wyear)
else
curterm = dble(chngpnt)
end if
go to 300
end if
if (curname.eq.'post1957 period' .or. curname.eq.
+ 'post1959 period' .or. curname.eq.'post1961 period') then
read (curname(7:8),'(i2)') chngpnt
if (wyear.gt.chngpnt) then
curterm = dble(wyear - chngpnt)
else
curterm = 0.0d0
end if
go to 300
end if
************************************************************************
* *
* This block of instructions is in here to retain compatibility with *
* earlier versions of this subroutine. *
* *
* We might wish to model the period effect by a surrogate variable. *
* Based on some exploration, we have chosen to model the period *
* effect by the following function: *
* *
* per59(wyear) = primary enrollment(wyear) - primary enrollment(59) *
* for wyear >= 59 and 0 otherwise, *
* *
* where 'primary enrollment' is defined to be the proportion of the *
* Iranian population attending primary education during the current *
* exposure-year. The primeduc vector contains these proportions. *
* *
************************************************************************
if (curname.eq.'per59') then
if (wyear.ge.59) then
curterm = primeduc(iyear) - primeduc(38)
else
curterm = 0.0d0
end if
go to 300
end if
************************************************************************
* *
* We might wish to model the period effect by a surrogate variable. *
* Based on some exploration, we have chosen to model the period *
* effect by functions of the following form (where yy = a two digit *
* year of an assumed changepoint, with 40 <= yy <= 77): *
* *
* post(wyear) = primary enrollment(wyear) - primary enrollment(yy)*
* for wyear >= yy and 0 otherwise, *
* *
* where 'primary enrollment' is defined to be the proportion of the *
* Iranian population attending primary education during the current *
* exposure-year. The primeduc vector contains these proportions. *
* *
************************************************************************
if (curname(1:6).eq.'post19' .and. curname(9:18).eq.
+ ' primary') then
twochar = curname(7:8)
if (twochar.ge.'40' .and. twochar.le.'77') then
read (twochar,'(i2)',err=240) chngpnt
if (wyear.ge.chngpnt) then
curterm = primeduc(iyear) - primeduc(chngpnt-21)
else
curterm = 0.0d0
end if
go to 300
end if
end if
240 continue
************************************************************************
* *
* If curname is the 'primary enrollment' identifier, include a *
* variable defined to be the proportion of the Iranian population *
* attending primary education during the current exposure-year. The *
* primeduc vector contains these proportions (which were determined *
* using data from B.R. Mitchell's "International Historical *
* Statistics: Africa and Asia"). *
* *
************************************************************************
if (curname.eq.'primary enrollment') then
curterm = primeduc(iyear)
go to 300
end if
************************************************************************
* *
* If curname is the 'primarybyprov' identifier, include a variable *
* defined to be the mean proportion of Iranian children attending *
* primary education for the province in which the woman resides (the *
* mean here is found by adding the proportion of female children to *
* the proportion of male children and dividing by 2). The data used *
* are in the femalebyprov and malebyprov vectors which contain values *
* from the 1976 Iranian census as provided by Akbar Aghajanian in *
* email dated May 24, 1992. *
* *
************************************************************************
if (curname.eq.'primarybyprov') then
read (provnumb,'(i2)') pvn
curterm = (femalebyprov(pvn) + malebyprov(pvn)) / 2.0d0
go to 300
end if
************************************************************************
* *
* If curname is the 'step19xx' identifier, include a variable taking *
* the value 1 for years 19xx and later and the value 0 for years *
* before 19xx, where xx is any number between 00 and 77, inclusive. *
* *
* As an example, the 'step1967' identifier could be used to include *
* an indicator variable which is 1 for woman exposure years beginning *
* with 1967 (the year in which the family protection act became law) *
* and later and 0 otherwise. *
* *
************************************************************************
if (curname(1:6).eq.'step19' .and. curname(9:18).eq.' ') then
twochar = curname(7:8)
if (twochar.ge.'00' .and. twochar.le.'77') then
read (twochar,'(i2)',err=260) chngpnt
curterm = compare('ge',wyear,chngpnt)
go to 300
end if
end if
260 continue
************************************************************************
* *
* If curname is the 'previous living' identifier, include previous in *
* the calculation of the model. Previous is a binary variable where *
* 1 -> previous child is still living at the time of the survey and *
* 0 -> previous child is no longer living at the time of the survey. *
* *
************************************************************************
if (curname.eq.'previous living') then
read (previous,'(f2.0)') curterm
go to 300
end if
************************************************************************
* *
* If curname is the 'previous this year' identifier, include a *
* covariate which equals previous except when previous is 0 (so the *
* previous child had died by time of survey) and prevlivy is not NA. *
* In the latter case, this covariate should be set to 1 only if the *
* previous child was still alive during the current exposure year. *
* *
************************************************************************
if (curname.eq.'previous this year') then
if (previous.eq.'00' .and. prevlivy.ne.'NA') then
read (prevlivy,'(i2)') yearsliv
curterm = compare('le',duration,yearsliv)
else
read (previous,'(f2.0)') curterm
end if
go to 300
end if
************************************************************************
* *
* If curname is any of the 'size' identifiers, include the corres- *
* ponding size covariate. Other than 'size' by itself, the five *
* size covariates available are identified by the names: *
* size0, size1, size2, size3 and size4 *
* *
* If curname is any of these five identifiers, include the corres- *
* ponding indicator variable in the calculation of the model. *
* *
************************************************************************
if (curname(1:4).eq.'size' .and. curname(6:18).eq.' ') then
read (size,'(i2)') isize
onechar = curname(5:5)
if (onechar.eq.' ') then
curterm = dble(isize)
go to 300
end if
if (onechar.ge.'0' .and. onechar.le.'4') then
read (onechar,'(i1)') digit
curterm = compare('eq',isize,digit)
go to 300
end if
end if
************************************************************************
* *
* If curname is the 'wedcat' identifier, include wedcat in the *
* calculation of the model. *
* *
************************************************************************
if (curname.eq.'wedcat') then
read (wedcat,'(f2.0)') curterm
go to 300
end if
************************************************************************
* *
* If curname is the 'hedcat' identifier, include hedcat in the *
* calculation of the model. *
* *
************************************************************************
if (curname.eq.'hedcat') then
read (hedcat,'(f2.0)') curterm
go to 300
end if
************************************************************************
* *
* The woman's age at marriage may be included in the model in either *
* of two fashions. Her marriage age may be included as a single *
* covariate. If curname is the 'agem' identifier, then include the *
* woman's marriage age in its original reported form. *
* *
************************************************************************
if (curname.eq.'agem') then
read (agem,'(f2.0)') curterm
go to 300
end if
************************************************************************
* *
* Instead of a single agem covariate, an alternative is to include *
* age at marriage as two covariates. The first of these is the *
* negative of the number of years before age 12 at which the woman *
* was married or 0 if the woman was married at age 12 or later. The *
* second of these is the number of years after age 12 at which the *
* woman was married or 0 if the woman was married before age 12. *
* *
* In curname these 2 covariates are referred to as 'early.agem' and *
* 'late.agem' (respectively). *
* *
************************************************************************
if (curname.eq.'early.agem') then
read (earlyage,'(f2.0)') curterm
go to 300
end if
if (curname.eq.'late.agem') then
read (lateage,'(f2.0)') curterm
go to 300
end if
************************************************************************
* *
* In question 104 of the individual survey, the woman was asked what *
* type of community she grew up in, up to age 12. Allowable answers *
* were: 1 -> city, or 2 -> village. This covariate is to be *
* included in the calculation if curname is the q104 identifier. *
* *
************************************************************************
if (curname.eq.'q104') then
read (q104,'(f2.0)') curterm
go to 300
end if
************************************************************************
* *
* There are currently ten available contraception covariates. Eight *
* of these are variables read in directly from covar.output (hence, *
* are copied right from the IFS database on the IBM VM machine). The *
* other two are variables which were constructed during the earlier *
* call to covarin. All ten variables are binary, where 0 -> woman *
* has not used this form of contraception and 1 -> woman has used *
* this form of contraception. *
* *
* The ten contraception covariates available are identified by the *
* names (the first eight are those from the IFS database): *
* pill, iud, other method, douche, condom, rhythm, withdrawal, *
* abstinence, modern, and traditional. *
* *
************************************************************************
if (curname.eq.'pill') then
read (pilluse,'(f2.0)') curterm
go to 300
end if
if (curname.eq.'iud') then
read (iud,'(f2.0)') curterm
go to 300
end if
if (curname.eq.'other method') then
read (otherm,'(f2.0)') curterm
go to 300
end if
if (curname.eq.'douche') then
read (douche,'(f2.0)') curterm
go to 300
end if
if (curname.eq.'condom') then
read (condom,'(f2.0)') curterm
go to 300
end if
if (curname.eq.'rhythm') then
read (rhythm,'(f2.0)') curterm
go to 300
end if
if (curname.eq.'withdrawal') then
read (withdraw,'(f2.0)') curterm
go to 300
end if
if (curname.eq.'abstinence') then
read (abstain,'(f2.0)') curterm
go to 300
end if
************************************************************************
* *
* Continue the contraception dummy covariate routines. *
* *
************************************************************************
if (curname.eq.'modern') then
read (modern,'(f2.0)') curterm
go to 300
end if
if (curname.eq.'traditional') then
read (tradition,'(f2.0)') curterm
go to 300
end if
************************************************************************
* *
* Since modern contraception was all but unheard of in Iran prior to *
* 1967, we might want to include a modern (or pill) contraception *
* covariate which has been scaled proportional to the overall use of *
* modern contraception in the population at large depending on the *
* precise exposure year. These two possible covariates are *
* identified by 'adjusted modern' and 'adjusted pill' in the covterms *
* matrix. *
* *
************************************************************************
if (curname.eq.'adjusted modern') then
curterm = usescale(modern,iyear)
go to 300
end if
if (curname.eq.'adjusted pill') then
curterm = usescale(pilluse,iyear)
go to 300
end if
************************************************************************
* *
* In regressions separate from this program we observed that the *
* relative proportion (relative to 1977) of contraception users in *
* Iran lagged 1 year was "roughly" related to the observed period *
* effect. If curname is the 'contraception lagx' identifier, *
* include this proportion of contraception users in Iran lagged by *
* x years, where x is any number between 0 and 9, inclusive. *
* *
************************************************************************
if (curname(1:17).eq.'contraception lag') then
onechar = curname(18:18)
if (onechar.ge.'0' .and. onechar.le.'9') then
read (onechar,'(i1)') digit
curterm = usescale('01',iyear - digit)
go to 300
end if
end if
************************************************************************
* *
* If curname is any of the 'province' identifiers, include the *
* corresponding province covariate. Other than 'province' by itself, *
* any of the 23 Iranian provinces, identified by numbers between 0 *
* 22 inclusive, may be included in the model using an indicator vari- *
* able set to a 1 if the woman associated with this row of intcovar *
* resided in the given province at the time of the survey or to a 0 *
* otherwise. Also, it is necessary that at least one province not be *
* specified in the set of covariates. All of the women who resided *
* in any of the non-included provinces together serve as a baseline *
* and as such are incorporated in the intercept. *
* *
* There are then a total of twenty-four possible province covariate *
* identifiers. They are province0, province1, ..., province9, *
* province10, ..., province19, province20, province21, province22, *
* and province. *
* *
************************************************************************
if (curname(1:8).eq.'province' .and. curname(11:18).eq.' ')
+ then
read (provnumb,'(i2)') pvn
if (curname(10:10).eq.' ') then
onechar = curname(9:9)
if (onechar.eq.' ') then
curterm = dble(pvn)
go to 300
end if
if (onechar.ge.'0' .and. onechar.le.'9') then
read (onechar,'(i1)') digit
curterm = compare('eq',pvn,digit)
go to 300
end if
else
twochar = curname(9:10)
if ( (twochar.ge.'10' .and. twochar.le.'19') .or.
+ (twochar.ge.'20' .and. twochar.le.'22') ) then
read (twochar,'(i2)') digit
curterm = compare('eq',pvn,digit)
go to 300
end if
end if
end if
************************************************************************
* *
* The ethnicity of individual women in the IFS is not available for *
* diplomatic reasons. However we do know the majority ethnicity for *
* each of the 23 Iranian provinces as determined during the 1976 *
* census. Akbar Aghajanian provided this information to us in email *
* dated May 24, 1992. There are 6 different majority ethnicities. *
* *
* If curname is any of the following, include an indicator variable *
* corresponding as to whether or not the selected ethnicity is the *
* majority ethnicity for the province in which the woman currently *
* resides: *
* arab, baluch, kurd, lor, persian or turk *
* *
************************************************************************
if (curname.eq.'arab') then
read (provnumb,'(i2)') pvn
curterm = compare('eq',ethnicity(pvn),0)
go to 300
end if
if (curname.eq.'baluch') then
read (provnumb,'(i2)') pvn
curterm = compare('eq',ethnicity(pvn),1)
go to 300
end if
if (curname.eq.'kurd') then
read (provnumb,'(i2)') pvn
curterm = compare('eq',ethnicity(pvn),2)
go to 300
end if
if (curname.eq.'lor') then
read (provnumb,'(i2)') pvn
curterm = compare('eq',ethnicity(pvn),3)
go to 300
end if
if (curname.eq.'persian') then
read (provnumb,'(i2)') pvn
curterm = compare('eq',ethnicity(pvn),4)
go to 300
end if
if (curname.eq.'turk') then
read (provnumb,'(i2)') pvn
curterm = compare('eq',ethnicity(pvn),5)
go to 300
end if
************************************************************************
* *
* There are currently sixteen available province level covariates, *
* all having been read in directly from province.output (hence, are *
* copied right from the IFS database on the IBM VM machine). All of *
* the province level covariates can be found in character form in the *
* provincb common block. They are generally decimal numbers with a *
* varying number of decimal places following the decimal point; some *
* are whole numbers. *
* *
* The sixteen province level covariates available are identified by *
* the names: *
* educ1, illratio, malelit, married, mlnonag, fmtoml, femsal, *
* electr, radio, tv, phone, urban3, uhincome, rhincome, clinirt *
* and userrt. *
* *
************************************************************************
if (curname.eq.'educ1') then
read (educ1,'(f4.1)') curterm
go to 300
end if
if (curname.eq.'illratio') then
read (illratio,'(f3.0)') curterm
go to 300
end if
if (curname.eq.'malelit') then
read (malelit,'(f4.1)') curterm
go to 300
end if
if (curname.eq.'married') then
read (married,'(f4.1)') curterm
go to 300
end if
if (curname.eq.'mlnonag') then
read (mlnonag,'(f4.1)') curterm
go to 300
end if
if (curname.eq.'fmtoml') then
read (fmtoml,'(f4.1)') curterm
go to 300
end if
if (curname.eq.'femsal') then
read (femsal,'(f4.1)') curterm
go to 300
end if
if (curname.eq.'electr') then
read (electr,'(f4.1)') curterm
go to 300
end if
if (curname.eq.'radio') then
read (radio,'(f4.1)') curterm
go to 300
end if
************************************************************************
* *
* Continue the province level covariate routines. *
* *
************************************************************************
if (curname.eq.'tv') then
read (tv,'(f4.1)') curterm
go to 300
end if
if (curname.eq.'phone') then
read (phone,'(f4.1)') curterm
go to 300
end if
if (curname.eq.'urban3') then
read (urban3,'(f4.1)') curterm
go to 300
end if
if (curname.eq.'uhincome') then
read (uhincome,'(f6.0)') curterm
curterm = curterm / 1.0d4
go to 300
end if
if (curname.eq.'rhincome') then
read (rhincome,'(f6.0)') curterm
curterm = curterm / 1.0d4
go to 300
end if
if (curname.eq.'clinirt') then
read (clinirt,'(f4.1)') curterm
curterm = curterm / 1.0d2
go to 300
end if
if (curname.eq.'userrt') then
read (userrt,'(f5.1)') curterm
curterm = curterm / 1.0d3
go to 300
end if
************************************************************************
* *
* This is the bottom of the covariate processing loop. It is an *
* error to fall through to this point without having recognized the *
* covariate identifier (that is, all of the covariate if-clauses end *
* with a goto to statement 300). So, if an error situation exists, *
* tell the user and set an error return code so our caller can know *
* what happened. *
* *
************************************************************************
write (0,*) 'unrecognizable covariate identifier encountered'
r15 = 4
return
300 product = product * curterm
400 continue
intcovar(icr,cn) = product
500 continue
************************************************************************
* *
* At this point the specified row of intcovar has been filled in. *
* *
************************************************************************
r15 = 0
return
end
************************************************************************
* *
* compare 4-16-92 *
* *
* This function determines whether or not the relation given by the *
* first argument holds between the second argument and the third *
* argument (the latter two must be integers). If the relationship *
* is true, this function returns a double precision value of 1.0d0. *
* If the relationship is false, the function returns a double *
* precision value of 0.0d0. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* operator = a character string containing one of the following *
* relational operators (they must be lower case): *
* 'eq', 'ne', 'lt', 'le', 'gt', or 'ge'. *
* *
* left = an integer containing the number which is to serve as *
* the lefthand argument to the relational operator. *
* *
* right = an integer containing the number which is to serve as *
* the righthand argument to the relational operator. *
* *
* *
* Outputs: *
* *
* compare = a double precision number which is set to 1.0d0 if the *
* relationship between the lefthand argument and the *
* righthand argument pertains. If the relationship *
* does not pertain, compare returns a value of 0.0d0. *
* If the operator is not one of the operators listed *
* above, then compare returns a value of -1.0d0. *
* *
************************************************************************
function compare(operator,left,right)
double precision compare
character operator*2
integer left
integer right
************************************************************************
* *
* For any valid form of the operator argument, perform the indicated *
* relational operation, and set compare to 1.0d0 if the relationship *
* is true, or 0.0d0 if the relationship is false. If the operator *
* argument is not one of the valid forms, set compare to -1.0d0. *
* *
************************************************************************
if ( operator .eq. 'eq' ) then
if (left.eq.right) then
compare = 1.0d0
else
compare = 0.0d0
end if
else if ( operator .eq. 'ne' ) then
if (left.ne.right) then
compare = 1.0d0
else
compare = 0.0d0
end if
else if ( operator .eq. 'lt' ) then
if (left.lt.right) then
compare = 1.0d0
else
compare = 0.0d0
end if
else if ( operator .eq. 'le' ) then
if (left.le.right) then
compare = 1.0d0
else
compare = 0.0d0
end if
else if ( operator .eq. 'gt' ) then
if (left.gt.right) then
compare = 1.0d0
else
compare = 0.0d0
end if
else if ( operator .eq. 'ge' ) then
if (left.ge.right) then
compare = 1.0d0
else
compare = 0.0d0
end if
else
compare = -1.0d0
end if
return
end
************************************************************************
* *
* dur3plus 9-23-92 *
* *
* Using the IFS dataset, we have derived a functional approximation *
* to the expected fertility rates for durations 3 and greater. The *
* following function returns this functional approximation given the *
* current value of the duration variable. The calculation of this *
* approximation is required at several different locations during the *
* calculations of the covariates to be stored in intcovar. In order *
* to facilitate any possible future changes in this approximation, I *
* chose to implement it as a function, to be called whenever it is *
* required. *
* *
* The constants used in this function were determined by performing *
* a weighted least squares regression fit of the natural logarithm of *
* expected fertility rate against duration, whose values were *
* empirically estimated from the IFS dataset. Because the covariate *
* matrix enters a logistic regression at the logit scale, it appears *
* most reasonable for this routine to return its approximation also *
* on the logit scale. *
* *
************************************************************************
function dur3plus(duration)
double precision dur3plus
integer duration
double precision rate
rate = exp(duration * -0.22934d0) * 0.72644d0
dur3plus = rate / (1.0d0 - rate)
return
end
************************************************************************
* *
* usescale 8-03-94 *
* *
* The usescale function multiplies whichever contraception-ever-used *
* variable is passed to it by the proportion of users of modern *
* contraception in the current exposure year and divides this product *
* by the proportion of users of modern contraception in the year of *
* the Iran Fertility Survey (using ratios determined from sources *
* independent of the IFS as reported in the 1991 paper by Akbar *
* Aghajanian titled "Social change, fertility and fertility change *
* in Iran"). *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* everuse = a character string containing one of the following *
* contraceptive ever used responses: *
* '00', '01', or 'NA'. *
* *
* year = an integer representing the current exposure year *
* (in the Iranian calendar, hence between 0 and 56, *
* inclusive). *
* *
* *
* Outputs: *
* *
* usescale = a double precision number between 0.0d0 and 1.0d0. *
* It is determined by taking the product of everuse *
* times the proportion of users of modern contraception *
* in the current exposure year divided by the proportion *
* of users of modern contraception in 1977 (1356 of *
* the Iranian calendar). If everuse was 'NA' or if *
* year was not between 0 and 56, then usescale returns *
* a -1.0d0. *
* *
************************************************************************
function usescale(everuse,year)
double precision usescale
character everuse*2
integer year
double precision scale(0:56) / 46*0.0d0, 0.009091d0, 0.254546d0,
+ 0.436364d0, 0.609091d0, 0.754546d0, 0.872727d0, 0.854546d0,
+ 0.845455d0, 0.836364d0, 0.945455d0, 1.0d0 /
double precision useinput
if ( ( everuse.eq.'00' .or. everuse.eq.'01' ) .and.
+ ( year.ge.0 .and. year.le.56 ) ) then
read (everuse,'(f2.0)') useinput
usescale = scale(year) * useinput
else
usescale = -1.0d0
end if
return
end
************************************************************************
* *
* stdindep 6-18-93 *
* *
* This subroutine performs a correlation transformation of an input *
* data matrix (a matrix in which the rows are observations of multi- *
* variate data and the columns represent variables). This subroutine *
* also returns the mean and standard deviation for each of the *
* variables present in the input data matrix. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* indep = a double precision matrix containing observed values *
* of all independent variables. Each column represents *
* a separate variable and each row is one complete *
* observation of all variables. This matrix should be *
* maxobs rows by numvars columns in dimension. The *
* actual number of rows used is passed as numobs. *
* *
* maxobs = an integer containing the allocated (i.e., maximum *
* available) number of rows in the input independent *
* variables matrix, indep. *
* *
* numvars = an integer containing the number of variables present *
* for each observation in the input independent *
* variables matrix, indep. In other words, numvars is *
* the number of columns of meaningful data in indep. *
* *
* numobs = an integer containing the number of observations (i.e. *
* number of rows) of data which has been passed in the *
* input independent variables matrix, indep. This *
* number must be less than or equal to maxobs. *
* *
* *
* Outputs: *
* *
* standard = a double precision matrix in which this subroutine is *
* to return the transformed version of the input matrix. *
* It must have the same dimensions as allocated for the *
* the input independent variables matrix, indep. This *
* subroutine will return the transformed data in exactly *
* the same order as the data is passed to it in the *
* input matrix. This subroutine is implemented in such *
* a fashion so that it will execute properly even if *
* the storage allocated for the input matrix coincides *
* with the storage for the output matrix. This permits *
* the calling routine to use one and the same matrix for *
* both input to and output from this subroutine so long *
* as the input values are not needed after calling this *
* subroutine. However proper execution of this subrou- *
* tine is no longer assured if the input matrix and the *
* output matrix only partially overlap each other. The *
* results will depend upon exactly how the two overlap. *
* *
************************************************************************
************************************************************************
* *
* Outputs (continued): *
* *
* mean = a double precision vector in which this subroutine is *
* to return the means of each of the variables passed in *
* the columns of the input independent variables matrix. *
* This vector must contain at least numvars elements. *
* *
* sd = a double precision vector in which this subroutine is *
* to return the standard deviations of each of the *
* variables passed in the columns of the input *
* independent variables matrix. This vector must *
* contain at least numvars elements. *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = numobs was less than 1 or greater than maxobs. *
* 8 = numvars was less than 1. *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine stdindep(indep,maxobs,numvars,numobs,standard,mean,sd,
+ r15)
double precision indep(maxobs,numvars)
integer maxobs
integer numvars
integer numobs
double precision standard(maxobs,numvars)
double precision mean(numvars)
double precision sd(numvars)
integer r15
************************************************************************
* *
* The following variables hold various temporary values used in this *
* subroutine. This includes any do-loop counters and similar such *
* temporary subscripts and indices. *
* *
* rootnm1 - square root of one less than number of observations *
* *
************************************************************************
double precision rootnm1
integer o1
integer v1
************************************************************************
* *
* Make sure that numobs is at least 1 and less than or equal to *
* maxobs. If not, return with error code 4. *
* *
************************************************************************
if (numobs.lt.1 .or. numobs.gt.maxobs) then
r15 = 4
return
end if
************************************************************************
* *
* Make sure that numvars is at least 1. If not, return with error *
* code 8. *
* *
************************************************************************
if (numvars.lt.1) then
r15 = 8
return
end if
************************************************************************
* *
* Find the sum of each column of indep and the sum of the squares of *
* each column of indep. The mean and sd vector arguments may be used *
* at least temporarily for this purpose. Of course it is important *
* to initialize these two vectors to zero before doing the addition. *
* *
************************************************************************
do 100 v1=1,numvars
mean(v1) = 0.0d0
sd(v1) = 0.0d0
100 continue
do 200 o1=1,numobs
do 150 v1=1,numvars
mean(v1) = indep(o1,v1) + mean(v1)
sd(v1) = indep(o1,v1) * indep(o1,v1) + sd(v1)
150 continue
200 continue
************************************************************************
* *
* Complete the calculation of the mean of each of the variables by *
* dividing the totals already in the mean vector by the number of *
* observations. *
* *
************************************************************************
do 300 v1=1,numvars
mean(v1) = mean(v1) / numobs
300 continue
************************************************************************
* *
* Next, calculate each variable's standard deviation multiplied by *
* the square root of one less than the number of observations. These *
* are the quantities which are actually needed during the calculation *
* of the correlation transformation which follows. *
* *
************************************************************************
do 400 v1=1,numvars
sd(v1) = dsqrt( sd(v1) - ( mean(v1) * mean(v1) * numobs ) )
400 continue
************************************************************************
* *
* Standardize each element of the input data matrix by subtracting *
* its mean and dividing by the quantities currently held in the sd *
* vector. Remember that these quantities are standard deviations *
* multiplied by the square root of one less than the number of *
* observations. *
* *
************************************************************************
do 600 o1=1,numobs
do 500 v1=1,numvars
standard(o1,v1) = ( indep(o1,v1) - mean(v1) ) / sd(v1)
500 continue
600 continue
************************************************************************
* *
* Finally, complete the calculation of the standard deviation of each *
* of the variables by dividing the values already in sd by the square *
* root of one less than the number of observations. *
* *
************************************************************************
rootnm1 = dsqrt( dble(numobs) - 1 )
do 700 v1=1,numvars
sd(v1) = sd(v1) / rootnm1
700 continue
************************************************************************
* *
* The correlation transformation has been completed successfully. So *
* return to the calling program with a return code indicating that no *
* problems were detected. *
* *
************************************************************************
r15 = 0
return
end
************************************************************************
* *
* wrtintc 2-14-94 *
* *
* This subroutine writes out the complete set of covariate data for *
* all woman exposure-years so long as there is no missing data for *
* the variables actually included in the model. The covariate data *
* is dumped to the disk file named in the dumpname variable (which is *
* located in the options common block). *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* intcovar = a double precision matrix containing the complete set *
* of covariate data for all women exposure-years in the *
* sample. This matrix should be maxint rows by numcovar *
* columns in dimension. The actual number of rows used *
* is passed as yearcnt. *
* *
* maxint = an integer containing the allocated (i.e., maximum *
* available) number of rows in the covariates matrix, *
* intcovar. *
* *
* numcovar = an integer containing the number of covariates passed *
* for each exposure-year (i.e., the number of columns in *
* the covariates matrix, intcovar. *
* *
* yearcnt = an integer containing the number of exposure-years of *
* covariate data which has been passed in the covariates *
* matrix, intcovar. This number must be less than or *
* equal to maxint. *
* *
* womenids = an integer vector containing the woman id associated *
* with each row of intcovar. Hence, it must contain *
* at least yearcnt number of elements. *
* *
* *
* Outputs: *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = yearcnt was less than 1 or greater than maxint. *
* 8 = numcovar was less than 1. *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine wrtintc(intcovar,maxint,numcovar,yearcnt,womenids,r15)
double precision intcovar(maxint,numcovar)
integer maxint
integer numcovar
integer yearcnt
integer womenids(yearcnt)
integer r15
************************************************************************
* *
* The fortran commands defining the option common block are in a *
* separate include file, 'optioncb.h'. Use of this include file is *
* not in accordance with the Fortran 77 Ansi standard. *
* *
************************************************************************
include 'optioncb.h'
************************************************************************
* *
* The following variables hold various temporary values used in this *
* subroutine. This includes any do-loop counters and similar such *
* temporary subscripts and indices. *
* *
* digits - number of digits past decimal point to show in output *
* *
* lincnt - number of lines to print for each exposure-year *
* valcnt - number of covariate values to include in each line *
* first - subscript of first covariate to dump to current line *
* last - subscript of last covariate to dump to current line *
* decr - decrement to first value in the line = (valcnt-1) *
* *
************************************************************************
character digits*1
integer lincnt
integer valcnt
integer first
integer last
integer decr
integer y1
integer l1
integer c1
************************************************************************
* *
* Make sure that yearcnt is at least 1 and less than or equal to *
* maxint. If not, return with error code 4. *
* *
************************************************************************
if (yearcnt.lt.1 .or. yearcnt.gt.maxint) then
r15 = 4
return
end if
************************************************************************
* *
* Make sure that numcovar is at least 1. If not, return with error *
* code 8. *
* *
************************************************************************
if (numcovar.lt.1) then
r15 = 8
return
end if
************************************************************************
* *
* Now dump the entire contents of the intcovar matrix to the file *
* named in dumpname (which is located in the options common block). *
* The first record of the dump file contains the total number of rows *
* and columns actually used in intcovar. The remaining records of *
* the dump file contain the covariates in row major order (i.e., the *
* first row precedes the second, and so on). The first row of each *
* woman's covariates begins with her woman id number; the id number *
* is not output in any of the subsequent rows for the same woman. *
* *
* The elements within each row of intcovar are dumped to the dump *
* file in as many records as are required, which depends on how many *
* significant digits past the decimal point are to be printed. This *
* number is controlled by the first character of the print options *
* flag (located in the options common block), which is copied to a *
* local variable, digits, for use in this subroutine. The number of *
* covariate values to dump to each record of the dump file (based on *
* the setting of the print options flag) is stored in valcnt. Hence, *
* each record of the dump file will contain valcnt double precision *
* numbers (i.e., valcnt columns of data, not counting the woman's id *
* number when it is present). Subsequent rows of intcovar start on *
* new output records. *
* *
************************************************************************
open (unit=41,file=dumpname,status='new')
write (41,'(i5,1x,i2)') yearcnt, numcovar
digits = printopt(1:1)
if (digits.eq.'3') then
valcnt = 8
else if (digits.eq.'2') then
valcnt = 9
else if (digits.eq.'4') then
valcnt = 7
else
valcnt = 6
end if
decr = valcnt - 1
lincnt = (numcovar + decr) / valcnt
************************************************************************
* *
* During the following loop, the contents of intcovar are actually *
* dumped to the dump file. *
* *
************************************************************************
do 700 y1=1,yearcnt
do 600 l1=1,lincnt
last = l1 * valcnt
first = last - decr
if (last.gt.numcovar) last = numcovar
if (first.eq.1) then
if (digits.eq.'6') then
write (41,'(i4.4,6(1x,f11.6))') womenids(y1),
+ (intcovar(y1,c1), c1=first,last)
else if (digits.eq.'3') then
write (41,'(i4.4,8(1x,f8.3))') womenids(y1),
+ (intcovar(y1,c1), c1=first,last)
else if (digits.eq.'2') then
write (41,'(i4.4,9(1x,f7.2))') womenids(y1),
+ (intcovar(y1,c1), c1=first,last)
else if (digits.eq.'4') then
write (41,'(i4.4,7(1x,f9.4))') womenids(y1),
+ (intcovar(y1,c1), c1=first,last)
else if (digits.eq.'5') then
write (41,'(i4.4,6(1x,f10.5))') womenids(y1),
+ (intcovar(y1,c1), c1=first,last)
end if
else
if (digits.eq.'6') then
write (41,'(4x,6(1x,f11.6))') (intcovar(y1,c1),
+ c1=first,last)
else if (digits.eq.'3') then
write (41,'(4x,8(1x,f8.3))') (intcovar(y1,c1),
+ c1=first,last)
else if (digits.eq.'2') then
write (41,'(4x,9(1x,f7.2))') (intcovar(y1,c1),
+ c1=first,last)
else if (digits.eq.'4') then
write (41,'(4x,7(1x,f9.4))') (intcovar(y1,c1),
+ c1=first,last)
else if (digits.eq.'5') then
write (41,'(4x,6(1x,f10.5))') (intcovar(y1,c1),
+ c1=first,last)
end if
end if
600 continue
700 continue
************************************************************************
* *
* At this point the entire intcovar matrix has been dumped, so close *
* the dump file and return successfully to our caller. *
* *
************************************************************************
close (unit=41)
r15 = 0
return
end
************************************************************************
* *
* logit 6-11-93 *
* *
* Subroutine to compute logistic regression. *
* *
************************************************************************
* input
* b initial parameter values
* y(n) number of positive responses from each binomial sample
* ynobs(n) number of bernoulli trials for each binomial sample
* this can be one
* x(n,ip) covariates, some column must be all 1's to include
* an intercept
* in actual row dimension of x
* n number of distinct binomial observations
* ip number of covariates
* mit maximum number of iterations
* output
* b parameter estimates
* ll2 -2 * the difference in the loglikelihood function for
* the full model minus the model with only an intercept
* term. Should not be interpreted when there is no
* intercept term in the model, but other differences
* could be computed in this case.
* var ip(ip+1)/2 vector contains the inverses hessian of
* the negative log likelihood. Only the upper triangle
* is stored. It is stacked by columns.
* ier = -2 ip exceeds 80, must increase internal work space
* = -1 poorly conditioned hessian matrix
* = 0 successful execution
* >mit failed to converge in mit step
* mit number of iterations required
subroutine logit(b,ll2,var,y,ynobs,x,in,n,ip,mit,loglike,ier)
double precision var(ip*(ip+1)/2), loglike
double precision y(n),ynobs(n),x(in,ip),b(ip),ll2
integer in,n,ip,mit,ier
* iter iteration number
* ich flag for stopping criteria
* w(80) work vector
* xsc(80) first partial derivatives of the negative log
* likelihood
* u(6480) contains cholesky decomposition of hessian of
* the loglikelihood function, limits ip to dimension
* 50 or less
double precision w(80),xsc(80),u(6480)
integer i,iter,ich
* check that ip does not exceed 80
if (ip.gt.80) then
ier = -2
go to 999
end if
iter = 1
* compute hessian and score function
900 call logfun(y,ynobs,x,b,var,xsc,in,n,ip)
* compute cholesky decomposition of the hessian and
* inverse of the hessian multiplied by the score
call chol(var,ip,(ip*(ip+1))/2,u,ier)
if (ier.gt.0) then
ier=-20+ier
go to 999
end if
call cholmu(u,xsc,ip,(ip*(ip+1)),w,xsc)
* compute update of b and check for convergence
ich=0
do 20 i=1,ip
if (abs(xsc(i)).gt.(.1e-9)*(abs(b(i))+.1)) ich=1
b(i)=b(i)-xsc(i)
20 continue
* check for convergence and failure to converge
if (ich.eq.0) then
ier=0
mit=iter
call logfun(y,ynobs,x,b,var,xsc,in,n,ip)
call chol(var,ip,(ip*(ip+1))/2,u,ier)
call cholin(u,ip,(ip*(ip+1))/2,w,var)
call logl2(y,ynobs,x,b,ll2,in,n,ip,loglike)
go to 999
end if
iter=iter+1
if (iter.gt.mit) then
ier=iter
call logfun(y,ynobs,x,b,var,xsc,in,n,ip)
call chol(var,ip,(ip*(ip+1))/2,u,ier)
call cholin(u,ip,(ip*(ip+1))/2,w,var)
call logl2(y,ynobs,x,b,ll2,in,n,ip,loglike)
go to 999
end if
* repeat iteration
go to 900
999 return
end
************************************************************************
* *
* logfun *
* *
* Subroutine to compute score function and hessian from the negative *
* log likelihood function of logistic regression. *
* *
************************************************************************
* input
* y(n) number of positive responses in each binomial sample
* y(nobs) number of bernoulli trials for each binomial sample
* x(n,ip) covariates, double precision array
* b(ip) current parameter values, double precision vector
* in actual row dimension of x
* n number of observations
* ip number of covariates
* output
* var ip(ip+1)/2 vector contains the hessian of the negative
* log likelihood
* xsc(ip) vector contains the score function
subroutine logfun(y,ynobs,x,b,var,xsc,in,n,ip)
double precision var(ip*(ip+1)/2),xsc(ip)
double precision x(in,ip),y(n),b(ip),ynobs(n)
integer in,ip,n
* v1 used to reference stacked storage in var
* i,j utility indices
* xprob prob that y=1 for given x,b
* ts(100) store covariates from an observation to reduce paging
* temp temporary storage of lone intermediate calculation
* used soley to shorten the length of the lines of code
double precision xprob,ts(100),temp
integer i,j,k,v1
* set var and xsc initial values to 0
v1=0
do 10 j=1,ip
xsc(j)=0.
do 20 k=1,j
v1=v1+1
var(v1)=0.
20 continue
10 continue
* for each observation
do 30 i=1,n
* compute the prob y=1; make copy of covariates
xprob=0.
do 40 j=1,ip
ts(j)=x(i,j)
xprob=xprob+ts(j)*b(j)
40 continue
xprob=exp(xprob)
xprob=xprob/(1.+xprob)
* compute and accumulate derivatives
v1=0
do 50 j=1,ip
if (ts(j).eq.0.) then
v1=v1+j
go to 50
end if
xsc(j)=xsc(j)+(ynobs(i)*ts(j)*xprob)-(y(i)*ts(j))
do 60 k=1,j
v1=v1+1
temp=ynobs(i)*ts(j)*ts(k)*xprob*(1.-xprob)
var(v1)=var(v1)+temp
60 continue
50 continue
30 continue
return
end
************************************************************************
* *
* *
************************************************************************
subroutine logl2(y,ynobs,x,b,ll2,in,n,ip,loglike)
double precision x(in,ip),y(n),b(ip),yprob,xprob,ynobs(n),
+ ll2,xntot,const,loglike
integer in,ip,n,i,j,k
* Compute -2 * the difference in the loglikelihood function for
* the full model minus the model with only an intercept term
ll2=0
yprob=0.
xntot=0.
const = 1.
* For each observation, accumulate contribution to the
* negative loglikelihood function
do 30 i=1,n
* Compute the prob y=1; and keep track of the proportionality
* constant for the log-likelihood calculation.
* The following was modified so that log only needs to be
* called once.
if ((y(i) .eq. ynobs(i)) .or. (y(i) .eq. 0)) then
go to 90
else
if ( 2*y(i) .le. ynobs(i) ) then
do 50 k = 1, y(i)
const = const * (ynobs(i) - k + 1.)/k
50 continue
else
do 80 k = 1, ynobs(i) - y(i)
const = const * (ynobs(i) - k + 1.)/k
80 continue
end if
end if
90 yprob=yprob+y(i)
xntot=xntot+ynobs(i)
xprob=0.
do 40 j=1,ip
xprob=xprob+x(i,j)*b(j)
40 continue
ll2=ll2+ynobs(i)*log(1.+exp(xprob))-y(i)*xprob
30 continue
loglike = log(const) - ll2
* compute difference with constant model and multiply by 2
xprob=yprob/xntot
xprob=log(xprob/(1-xprob))
ll2=(xntot*log(1.+exp(xprob))-yprob*xprob)-ll2
ll2=2.*ll2
return
end
************************************************************************
* *
* cholesky routines *
* *
* The first routine computes the cholesky decompostion. *
* *
* Calculate cholesky decomposition of a and place result in U. *
* *
* Algorithm as 6 Applied Statistics, (1968) vol. 17, p.195. *
* *
************************************************************************
* a double precision array length nn that contains a positive
* definite matrix stored in symmetric storage mode with
* upper triangle stored by columns
* n dimension of a
* nn length of symmetrically stored matrices a and u
* u on output, contains cholesky decomposition of a
* ifault output:
* 1 if n is less than 1
* 2 if a is not positive definite
* 3 if nn is not equal to n(n+1)/2
* 0 otherwise
subroutine chol(a,n,nn,u,ifault)
double precision a(nn),u(nn),eta,eta2,x,w,zero,zabs,zsqrt
integer n,nn,ifault,j,k,l,m,kk,irow,icol,ii,i
* eta is used to determine if a is positive definite.
* for 32 bit representations, 1.0e-5 is appropriate,
* for 64 bit representations, 1.0e-9 can be used
data eta,zero /1.0e-5,0.0/
zabs(x)=abs(x)
zsqrt(x)=sqrt(x)
ifault=1
if (n.le.0) return
ifault=3
if (nn.ne.(n*(n+1)/2)) return
ifault=2
j=1
k=0
eta2=eta*eta
ii=0
do 80 icol=1,n
ii=ii+icol
x=eta2*a(ii)
l=0
kk=0
do 40 irow=1,icol
kk=kk+irow
k=k+1
w=a(k)
m=j
do 10 i=1,irow
l=l+1
if (i.eq.irow) go to 20
w=w-u(l)*u(m)
m=m+1
10 continue
20 if (irow .eq. icol) go to 50
if (u(l) .eq. zero) go to 30
u(k)=w/u(l)
go to 40
30 if (w*w .gt. zabs(x*a(kk))) return
u(k)=zero
40 continue
* modified to quit if A is not strictly positive definite
50 if (zabs(w) .le. zabs(eta*a(k))) return
if (w.lt.zero) return
u(k)=zsqrt(w)
j=j+icol
80 continue
ifault=0
return
end
************************************************************************
* *
* cholmu *
* *
* Compute A inverse * b where A is positive definite and its cholesky *
* decomposition is in U--- the main work here is in the forward and *
* backsolve algorithms. *
* *
************************************************************************
* u contains upper triangular cholesky decomposition of A
* b contains vector to be multiplied by the inverse
* n dimension of space--- length of b
* nn length of u=n(n+1)/2
* x product of inv(A) and b
* ga work vector that hold intermediate output
subroutine cholmu(u,b,n,nn,ga,x)
double precision u(nn),b(n),x(n),ga(n),sub
integer n,nn,i,j,ind,irow,ierc
* first forwardsolve the transpose of U with b. tran(U) will be
* lower triangular. Since u is stored like a symmetric matrix,
* using the sssub function applied to the lower half of U will
* return the transpose elements of U.
do 10 i=1,n
* work on rows from top to bottom
* accumulate products from col 1 to (i-1)
sub=0.
do 20 j=1,(i-1)
call sssub(n,i,j,ind,ierc)
sub=sub+u(ind)*ga(j)
20 continue
call sssub(n,i,i,ind,ierc)
ga(i)=(b(i)-sub)/(u(ind))
10 continue
* backsolve this intermediate solution to obtain final solution
do 30 i=1,n
* work on rows from bottom to top
irow=(n+1)-i
* accumulate products for (irow+1) to p elements
sub=0
do 40 j=(irow+1),n
call sssub(n,irow,j,ind,ierc)
sub=sub+u(ind)*x(j)
40 continue
call sssub(n,irow,irow,ind,ierc)
x(irow)=(ga(irow)-sub)/u(ind)
30 continue
return
end
************************************************************************
* *
* cholin *
* *
* Compute A inverse where A is positive definite and its cholesky *
* decomposition is in U. The routine does not check that the matrix *
* is invertible. This should be done using the output of the *
* cholesky routine. *
* *
* This is a modification of Algorithm 7, Applied Statistics (1968) *
* vol. 17, p. 198. *
* *
************************************************************************
* u contains upper triangular cholesky decomposition of A
* n dimension of space--- length of b
* nn length of u=n(n+1)/2
* w work vector
* c inv(A)
subroutine cholin(u,n,nn,w,c)
double precision u(nn),c(nn),w(n),x,zero,one
integer n,nn,i,k,l,irow,ndiag,icol,jcol,mdiag
data zero,one /0.0,1.0/
* place u in c for inplace inversion to remain consistent
* with the algorithm
do 5 i=1,nn
c(i)=u(i)
5 continue
irow=n
ndiag=nn
10 l=ndiag
do 20 i=irow,n
w(i)=c(l)
l=l+i
20 continue
icol=n
jcol=nn
mdiag=nn
30 l=jcol
x=zero
if (icol.eq.irow) x=one/w(irow)
k=n
40 if (k.eq.irow) go to 50
x=x-w(k)*c(l)
k=k-1
l=l-1
if (l.gt.mdiag) l=l-k+1
go to 40
50 c(l)=x/w(irow)
if (icol.eq.irow) go to 80
mdiag=mdiag-icol
icol=icol-1
jcol=jcol-1
go to 30
80 ndiag=ndiag-irow
irow=irow-1
if (irow.ne.0) go to 10
return
end
************************************************************************
* *
* sssub 6-06-93 *
* *
* This subroutine converts symmetric square matrix references to *
* references for symmetric storage mode. Lower triangular references *
* are mapped into upper triangular references. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* dimen = an integer containing the dimension (i.e. the number *
* of rows and number of columns) of the square matrix. *
* *
* rownum = an integer containing the row subscript of the element *
* being referenced. *
* *
* colnum = an integer containing the column subscript of the *
* element being referenced. *
* *
* *
* Outputs: *
* *
* result = an integer which will be set to the single subscript *
* into a symmetric storage mode vector which is the *
* equivalent of the element subscripted by the input row *
* and column subscripts. *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = rownum was less than 1 or greater than dimen. *
* 8 = colnum was less than 1 or greater than dimen. *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine sssub(dimen,rownum,colnum,result,r15)
integer dimen
integer rownum
integer colnum
integer result
integer r15
************************************************************************
* *
* Make sure that rownum is at least 1 and less than or equal to *
* dimen. If not, return with error code 4. *
* *
************************************************************************
if (rownum.lt.1 .or. rownum.gt.dimen) then
r15 = 4
return
end if
************************************************************************
* *
* Make sure that colnum is at least 1 and less than or equal to *
* dimen. If not, return with error code 8. *
* *
************************************************************************
if (colnum.lt.1 .or. colnum.gt.dimen) then
r15 = 8
return
end if
************************************************************************
* *
* Now calculate the symmetric storage mode subscript equivalent for *
* the square matrix element referred to by the rownum and colnum *
* arguments to this subroutine. *
* *
************************************************************************
if (rownum.gt.colnum) then
result = colnum + ( rownum - 1 ) * rownum / 2
else
result = rownum + ( colnum - 1 ) * colnum / 2
end if
r15 = 0
return
end
************************************************************************
* *
* wrtfixed 10-24-94 *
* *
* This subroutine writes the final results of the logistic regression *
* to the disk file named in the logitname variable (which is located *
* in the options common block). The output contains a complete set *
* of estimated fixed effect parameters with their standard errors and *
* t-values followed by model information such as the parity selection *
* criteria, the log-likelihood of the model, the corresponding value *
* of Bic, the total number of women exposure-years, the total number *
* of births and the error return code from the logit subroutine. *
* These are optionally followed by a complete output of the estimated *
* correlation matrix of the fixed effect parameters. *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* varmat = a double precision square matrix which contains the *
* estimated covariances of the fixed effect parameter *
* estimates. The number of rows and columns in varmat *
* is passed in the maxcovar argument. *
* *
* maxcovar = an integer containing the maximum number of covariates *
* which may be included in the model. The covariance *
* matrix of the parameter estimates, varmat, contains *
* this many rows and columns. *
* *
* fixed = a double precision vector containing the final *
* estimates of the fixed effect parameters. There are *
* numcovar elements in this vector. *
* *
* numcovar = an integer containing the actual number of covariates *
* included in the model. There are numcovar elements *
* in the fixed effect parameters vector, fixed. *
* *
* covnames = a vector of character strings with each string con- *
* taining a covariate identifier in precisely the format *
* read from the model input file. These identifiers *
* are included in the output produced by wrtfixed. *
* There are numcovar covariate identifiers in covnames. *
* *
* parint = an integer containing the first (chronologically *
* earliest) parity included in the model. *
* *
* partype = a one character character string set to either a blank *
* or a '+' character depending on which is read from the *
* model input file. *
* *
* ll2 = a double precision number which contains -2 * the *
* difference in the log-likelihood function for the full *
* model minus the model with only an intercept term. *
* *
* llk = a double precision number which contains the log- *
* likelihood of the full model (i.e., the model with *
* the full complement of parameters included). *
* *
************************************************************************
************************************************************************
* *
* Inputs (continued): *
* *
* yearcnt = an integer containing the actual total number of *
* women exposure-years used in estimating the model. *
* *
* birthcnt = an integer containing the total number of births to *
* all women in the sample. *
* *
* itercnt = an integer containing the actual number of iterations *
* it took for the logit subroutine to converge to the *
* final estimates. *
* *
* logitrc = an integer containing the error return code from the *
* logit subroutine. *
* *
* *
* Outputs: *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = numcovar was less than 1 or greater than 80. *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine wrtfixed(varmat,maxcovar,fixed,numcovar,covnames,
+ parint,partype,ll2,llk,yearcnt,birthcnt,itercnt,logitrc,r15)
double precision varmat(maxcovar,maxcovar)
integer maxcovar
double precision fixed(numcovar)
integer numcovar
character covnames(numcovar)*(*)
integer parint
character partype*1
double precision ll2
double precision llk
integer yearcnt
integer birthcnt
integer itercnt
integer logitrc
integer r15
************************************************************************
* *
* The fortran commands defining the option common block are in a *
* separate include file, 'optioncb.h'. Use of this include file is *
* not in accordance with the Fortran 77 Ansi standard. *
* *
************************************************************************
include 'optioncb.h'
************************************************************************
* *
* The following variables hold various temporary values used in this *
* subroutine. This includes any do-loop counters and similar such *
* temporary subscripts and indices. *
* *
* correlat - estimated correlations of the fixed effect parameters *
* (correlat must contain at least numcovar elements) *
* stddev - standard deviation of the current parameter estimate *
* tstat - estimated t-value of the current parameter estimate *
* bic - this model's Bic (Bayes information criteria) score *
* *
* digits - number of digits past decimal point to show in output *
* *
* lastsig - index to the last non-blank character of provset *
* lnblnk - Sun subroutine which returns the index of the last *
* non-blank character *
* *
************************************************************************
double precision correlat(80)
double precision stddev
double precision tstat
double precision bic
character digits*1
integer lastsig
integer lnblnk
integer r1
integer c1
************************************************************************
* *
* Make sure that numcovar is at least 1 and less than or equal to 80 *
* so that we don't attempt to write past the end of the correlat *
* vector. If not, return with error code 4. *
* *
************************************************************************
if (numcovar.lt.1 .or. numcovar.gt.80) then
r15 = 4
return
end if
************************************************************************
* *
* Open the logistic regression output file. Begin the output with *
* a couple lines of user specified parity and province selection *
* information. *
* *
************************************************************************
open (unit=42,file=logitname,status='new')
write (42,'(''Parity Selection = '',i1,a)') parint, partype
lastsig = lnblnk(provset)
write (42,'(''Province Selection = '',a)') provset(1:lastsig)
************************************************************************
* *
* Write the final parameter estimates, standard errors and t values. *
* The number of significant digits past the decimal point to be *
* printed is controlled by the first character of the print options *
* flag (located in the options common block), which is copied to a *
* local variable, digits, for use in this subroutine. *
* *
************************************************************************
write (42,'(/3x,''Parameter'',4x,''Covariate'',13x,''Estimate'',
+ 6x,''Std. Err'',6x,''t value'')')
digits = printopt(1:1)
do 300 c1=1,numcovar
stddev = dsqrt(varmat(c1,c1))
tstat = fixed(c1) / stddev
if (digits.eq.'6') then
write (42,'(4x,''b('',i2,'')'',5x,a,1x,f10.6,2(2x,f12.6))')
+ c1,covnames(c1)(1:21), fixed(c1), stddev, tstat
else if (digits.eq.'3') then
write (42,'(4x,''b('',i2,'')'',5x,a,1x,f7.3,2(5x,f9.3))')
+ c1,covnames(c1)(1:24), fixed(c1), stddev, tstat
else if (digits.eq.'2') then
write (42,'(4x,''b('',i2,'')'',5x,a,1x,f6.2,2(6x,f8.2))')
+ c1,covnames(c1)(1:25), fixed(c1), stddev, tstat
else if (digits.eq.'4') then
write (42,'(4x,''b('',i2,'')'',5x,a,1x,f8.4,2(4x,f10.4))')
+ c1,covnames(c1)(1:23), fixed(c1), stddev, tstat
else if (digits.eq.'5') then
write (42,'(4x,''b('',i2,'')'',5x,a,1x,f9.5,2(3x,f11.5))')
+ c1,covnames(c1)(1:22), fixed(c1), stddev, tstat
end if
300 continue
************************************************************************
* *
* Also print the following global statistics, including a calculation *
* of bic. *
* *
************************************************************************
bic = dlog( dble(yearcnt) ) * dble( numcovar - 1 ) - ll2
write (42,'(/''Likelihood Ratio Test Against Null Model = '',
+ f10.3)') ll2
write (42,'(''Log-Likelihood of the Model = '',f10.3)') llk
write (42,'(''Total Woman-Years of Exposure = '',i6)') yearcnt
write (42,'(''Total Number of Births = '',i6)') birthcnt
write (42,'(''Bic for this Model = '',f10.3)') bic
write (42,'(''Number of iterations to convergence = '',i2)')
+ itercnt
write (42,'(''Error = '',i4)') logitrc
************************************************************************
* *
* If the second character of the print options flag is a 'c', we are *
* to also print out the entire estimated correlation matrix of the *
* parameters. *
* *
************************************************************************
if (printopt(2:2).eq.'c') then
write(42,'(//''Estimated correlation matrix of the parameters:''
+ )')
do 700 r1=1,numcovar
do 600 c1=1,r1
correlat(c1) = varmat(r1,c1) / dsqrt( varmat(r1,r1) *
+ varmat(c1,c1) )
600 continue
write (42,'(40(1x,f5.2))') (correlat(c1), c1=1,r1)
700 continue
end if
close (unit=42)
r15 = 0
return
end
************************************************************************
* *
* wrtresid 8-03-94 *
* *
* This subroutine writes out the observed and fitted values along *
* with residuals for each woman exposure-year in the sample to the *
* disk file named in the residname variable (which is located in the *
* options common block). *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* intcovar = a double precision matrix containing the complete set *
* of covariate data for all women exposure-years in the *
* sample. This matrix should be maxint rows by numcovar *
* columns in dimension. The actual number of rows used *
* is passed as yearcnt. *
* *
* maxint = an integer containing the allocated (i.e., maximum *
* available) number of rows in the covariates matrix, *
* intcovar. *
* *
* numcovar = an integer containing the number of covariates passed *
* for each exposure-year. This must then also be the *
* number of elements in the fixed effects parameter *
* vector, fixed. *
* *
* yearcnt = an integer containing the number of exposure-years of *
* data to use in the calculation of the residuals. *
* This number must be less than or equal to maxint. *
* *
* fixed = a double precision vector containing the final *
* estimates of the fixed effect parameters. This vector *
* contains numcovar elements. *
* *
* random = a double precision vector containing the final *
* estimates of the woman specific random effect *
* parameters. This vector contains womaxcnt elements. *
* *
* womaxcnt = an integer containing the maximum number of women in *
* the sample. In other words, womaxcnt contains the *
* largest woman identification number used. It is the *
* number of elements in the woman specific random *
* effects vector, random. *
* *
* womenids = an integer vector containing the woman id associated *
* with each row of intcovar. Hence, it must contain *
* at least yearcnt number of elements. This vector is *
* used to subscript into the woman specific random *
* effects vector, random. *
* *
* response = a double precision vector with yearcnt elements where *
* each element is set to 1 if this woman gave birth *
* during the corresponding exposure-year and is set to *
* 0 if this woman did not give birth during the year. *
* Each element of the response vector is associated with *
* the corresponding row of the intcovar matrix. *
* *
************************************************************************
************************************************************************
* *
* Outputs: *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = yearcnt was less than 1 or greater than maxint. *
* 8 = numcovar was less than 1. *
* 12 = the response vector passed to wrtresid was *
* invalid (it must contain only 1's and 0's). *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine wrtresid(intcovar,maxint,numcovar,yearcnt,fixed,random,
+ womaxcnt,womenids,response,r15)
double precision intcovar(maxint,numcovar)
integer maxint
integer numcovar
integer yearcnt
double precision fixed(numcovar)
integer womaxcnt
double precision random(womaxcnt)
integer womenids(yearcnt)
double precision response(yearcnt)
integer r15
************************************************************************
* *
* The fortran commands defining the option common block are in a *
* separate include file, 'optioncb.h'. Use of this include file is *
* not in accordance with the Fortran 77 Ansi standard. *
* *
************************************************************************
include 'optioncb.h'
************************************************************************
* *
* The following variables hold various temporary values used in this *
* subroutine. This includes any do-loop counters and similar such *
* temporary subscripts and indices. *
* *
* exponent - current exposure-year inner product summation variable *
* residual - difference between the observed and fitted response *
* raised - exponentiated form of current inner product summation *
* fitted - model fitted value for the current exposure-year *
* *
************************************************************************
double precision exponent
double precision residual
double precision raised
double precision fitted
integer wid
integer y1
integer c1
************************************************************************
* *
* Make sure that yearcnt is at least 1 and less than or equal to *
* maxint. If not, return with error code 4. *
* *
************************************************************************
if (yearcnt.lt.1 .or. yearcnt.gt.maxint) then
r15 = 4
return
end if
************************************************************************
* *
* Make sure that numcovar is at least 1. If not, return with error *
* code 8. *
* *
************************************************************************
if (numcovar.lt.1) then
r15 = 8
return
end if
************************************************************************
* *
* Check that the response vector passed to this subroutine contains *
* only 1's and 0's. If it contains any other values, return to the *
* caller with an error code of 12. *
* *
************************************************************************
do 100 y1=1,yearcnt
if (response(y1).ne.1.0d0 .and. response(y1).ne.0.0d0) then
r15 = 12
return
end if
100 continue
************************************************************************
* *
* Now write the woman's id number, observed value, fitted value and *
* residual for each exposure-year to the file named in residname *
* (which is located in the options common block). *
* *
************************************************************************
open (unit=45,file=residname,status='new')
do 700 y1=1,yearcnt
wid = womenids(y1)
exponent = random(wid)
do 600 c1=1,numcovar
exponent = exponent + intcovar(y1,c1) * fixed(c1)
600 continue
raised = dexp(exponent)
fitted = raised / ( raised + 1.0d0 )
residual = response(y1) - fitted
write (45,'(i4.4,1x,i1,1x,f9.6,1x,f8.6,1x,f9.6)') wid,
+ int(response(y1)), exponent, fitted, residual
700 continue
close (unit=45)
r15 = 0
return
end
************************************************************************
* *
* handler routines 2-15-94 *
* *
* The subroutines included in this file are ieee exception handling *
* interface routines. They are grouped together here so that it will *
* be easy to convert them to no-operation subroutines in order that *
* any of the IFS programs can be run on machines which don't support *
* ieee floating point arithmetic. This can be done by commenting out *
* all the executable commands contained in the handlon and handloff *
* routines and inserting a command setting r15 to 0 in each of the *
* subroutines. The three subroutines included in this file are: *
* *
* handlon - called to begin trapping exceptions *
* handloff - called to stop trapping exceptions *
* handler - the actual exception handler itself *
* *
************************************************************************
************************************************************************
* *
* handlon *
* *
* Call handlon to begin trapping ieee floating point arithmetic *
* exceptions. This will usually be the first thing done in any main *
* program. The option argument is used to specify which type(s) of *
* exceptions should be trapped. There are five different types of *
* ieee floating point exceptions. These are 'inexact', 'division by *
* zero', 'underflow', 'invalid' and 'overflow'. Inexact exceptions *
* are usually of no consequence and so the caller will often want to *
* ignore these. In addition there are situations where underflow *
* exceptions are unavoidable and so once again the caller may want to *
* suppress these. The preceding explains the rational behind the *
* current implentation of three distinct exception trapping options. *
* The three options are *
* *
* 'all' = trap all exceptions *
* 'all-inexact' = trap all exceptions except inexact exceptions *
* 'common' = trap 'division by zero', 'invalid' and *
* 'overflow' exceptions only *
* *
* *
* Inputs: *
* *
* option = a character string containing one of the following *
* handling options: *
* 'all', 'all-inexact', or 'common' *
* *
* *
* Outputs: *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* 4 = unrecognizable option argument. *
* 8 = ieee.handler returned with a nonzero return *
* code. *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine handlon(option,r15)
character option*(*)
integer r15
************************************************************************
* *
* The following variables hold various constants and temporary values *
* used in this subroutine such as loop counters and so forth. *
* *
************************************************************************
external handler
integer ieee_handler
integer rc
************************************************************************
* *
* Check the handling option to determine which exceptions are to be *
* trapped. All exceptions will be trapped if the option is 'all'. *
* All exceptions except floating point inexact exceptions will be *
* trapped if the option is 'all-inexact'. Only 'division by zero', *
* 'invalid' and 'overflow' exceptions will be trapped if the option *
* is 'common'. If the option passed is not one of these options, *
* return to the caller with an error code of 4 (which is the default *
* setting for r15). *
* *
************************************************************************
r15 = 4
if (option(1:3).eq.'all') then
rc = ieee_handler('set', 'all', handler)
if (rc.ne.0) then
r15 = 8
return
end if
r15 = 0
end if
if (option.eq.'all-inexact') then
rc = ieee_handler('clear', 'inexact')
if (rc.ne.0) then
r15 = 8
return
end if
end if
if (option.eq.'common') then
rc = ieee_handler('set', 'common', handler)
if (rc.ne.0) then
r15 = 8
return
end if
r15 = 0
end if
return
end
************************************************************************
* *
* handloff *
* *
* Call handloff to terminate all trapping of ieee exceptions. This *
* will usually be the last thing done before stopping the main *
* program. The current implementation of handloff does not support *
* any options other than clearing all exception traps. *
* *
* *
* Inputs: *
* *
* option = a character string containing one of the following *
* handling options: *
* 'all' *
* *
* *
* Outputs: *
* *
* r15 = an integer valued error return code. This variable *
* is set to 0 if no errors were encountered. *
* Otherwise, r15 can assume the following values: *
* *
* -4 = unrecognizable option argument. *
* = ieee.flags returned with a nonzero return flag *
* of rc (see man ieee.flags for more details). *
* No other possible values are currently in use. *
* *
************************************************************************
subroutine handloff(option,r15)
character option*(*)
integer r15
************************************************************************
* *
* The following variables hold various constants and temporary values *
* used in this subroutine such as loop counters and so forth. *
* *
************************************************************************
character reply*12
external handler
integer ieee_handler
integer ieee_flags
integer rc
************************************************************************
* *
* Check if the handling option is valid. Currently only the 'all' *
* option is valid, so if it is anything else, return with an error *
* code of -4 (which is the default setting for r15). If the option *
* is 'all', then clear all ieee floating point exception handling. *
* *
************************************************************************
r15 = -4
if (option.eq.'all') then
r15 = ieee_flags('get', 'exception', 'all', reply)
rc = ieee_flags('clear', 'exception', 'all', reply)
rc = ieee_handler('clear', 'all')
end if
return
end
************************************************************************
* *
* handler *
* *
* This function is a general purpose user-specified signal handler. *
* When properly enabled by the user (see below), this function will *
* write a warning message to the standard error file whenever any *
* enabled ieee floating point arithmetic exception occurs. The *
* warning message contains the exception code number, a standard name *
* for the exception and the decimal address at which the exception *
* occurred. *
* *
* In order to enable this exception signal handler, the user should *
* include a sequence of commands such as the following as the first *
* executable commands in his/her fortran program. For more details *
* on the use of this exception signal handler see the chapter on *
* floating point arithmetic in the Sun Fortran User's Guide. *
* *
* Assuming that the user has defined rc to be an integer variable, *
* the following sequence of commands could be used to enable this *
* exception signal handler function. *
* *
* external handler *
* rc = ieee_handler('set', 'all', handler) *
* if (rc.ne.0) *
* *
************************************************************************
************************************************************************
* *
* Inputs: *
* *
* signumb = an integer containing the signal number for a floating *
* point exception. *
* *
* excpcode = an integer containing the address of the floating *
* point exception block for the type of exception which *
* has occurred. *
* *
* context = a five element integer vector containing data *
* pertinent to this exception. The fourth element *
* contains the address (current program counter value) *
* at which the exception occurred. *
* *
* *
* Outputs: < none > *
* *
************************************************************************
integer function handler(signumb,excpcode,context)
integer signumb
integer excpcode
integer context(5)
************************************************************************
* *
* The following variables hold various constants and temporary values *
* used in this subroutine such as loop counters and so forth. *
* *
************************************************************************
character excpname(5)*16 / 'inexact', 'division by zero',
+ 'underflow', 'invalid', 'overflow' /
integer excploc
integer loc
integer e1
************************************************************************
* *
* Fill in and print the floating point exception warning message. *
* *
************************************************************************
excploc = loc(excpcode)
e1 = (excploc - 192) / 4
write (0,'(''ieee exception code '',i3,'' ('',a,
+ '') occurred at pc '',i6)') excploc, excpname(e1), context(4)
return
end
SHAR_EOF
chmod 644 'eha.f'
fi
if test -f 'intcovcb.h'
then
echo shar: "will not over-write existing file 'intcovcb.h'"
else
cat << \SHAR_EOF > 'intcovcb.h'
************************************************************************
* *
* intcovcb 3-17-94 *
* *
* The large internal event history covariates matrix, intcovar, has *
* been located in its own uninitialized common block, intcovcb, in *
* order to not waste disk space for it in any executable created for *
* use with the Iran Fertility Study data. *
* *
************************************************************************
************************************************************************
* *
* The following five parameters are defined here so that the maximum *
* number of covariates supported in the implementation may be readily *
* altered, simply by modifying this include file. *
* *
************************************************************************
integer ppmaxcov
parameter (ppmaxcov=64)
integer ppmaxcp1
parameter (ppmaxcp1=ppmaxcov+1)
integer ppmaxcp3
parameter (ppmaxcp3=ppmaxcov+3)
integer ppmaxcp7
parameter (ppmaxcp7=ppmaxcov+7)
integer ppmatmax
parameter (ppmatmax=ppmaxcov*ppmaxcp1/2)
************************************************************************
* *
* Intcovar is preceded by two integer variables which will be set to *
* the number of rows and number of columns available in intcovar. *
* *
* maxint - allocated number of rows of the covariates matrix *
* maxcovar - the maximum number of covariates which can fit *
* intcovar - the matrix of discrete time unit covariates *
* *
************************************************************************
common /intcovcb/ maxint, maxcovar, intcovar
integer maxint
integer maxcovar
double precision intcovar(93100,ppmaxcov)
SHAR_EOF
chmod 644 'intcovcb.h'
fi
if test -f 'optioncb.h'
then
echo shar: "will not over-write existing file 'optioncb.h'"
else
cat << \SHAR_EOF > 'optioncb.h'
************************************************************************
* *
* optioncb 10-28-94 *
* *
* The options common block contains a number of globally referenced *
* variables, which are either character strings or integers. The *
* character strings are mostly either external filenames (up to 24 *
* characters long) or contain 'T'rue or 'F'alse indicators of which *
* global option(s) the user has selected. The integers are generally *
* information on number of iterations, seed values for the random *
* number generation subroutines, or other such similar counts. *
* *
************************************************************************
common /optioncb/ optinname, modelname, fixdname, randname,
+ dumpname, logitname, gibbsname, residname, mlfrprioname,
+ printopt, provset, fixedmethod, fixedprior, fixedsupport,
+ randprior, varprior, fixedmargvar, varshape, varscale,
+ firstvariance, iteracnt, burnin, priosize, parmdps, cseed,
+ tseed, randtracecnt, fixedcnt, randcnt, fixedndx, womanid,
+ ieee, standardize, intcdump, residout, gibbstrace, gibbstime,
+ mlfrpriotrace, fixedchk, randchk, varchk, keep4707, longeropen
************************************************************************
* *
* Optinname is the name of the options input file which was passed *
* from the command line as an argument to the wfsfit or wfsgibbs *
* programs. It is the name of the file which contains the user's *
* settings for the other options in this common block. The default *
* name is ' '. *
* *
************************************************************************
character optinname*24
************************************************************************
* *
* Modelname is the name of the input file from which the list of *
* variables making up the current model is read. The default name is *
* 'ifsmodel'. *
* *
************************************************************************
character modelname*24
************************************************************************
* *
* Fixdname is the name of the input file containing the initial *
* values to set the fixed effect parameters to. The default name is *
* 'fixdinit'. *
* *
************************************************************************
character fixdname*24
************************************************************************
* *
* Randname is the name of the input file containing the initial *
* values to set the random effect parameters to. The default name is *
* 'randinit'. *
* *
************************************************************************
character randname*24
************************************************************************
* *
* Dumpname is the name to use for the output file created when the *
* intcdump option is selected. The default name is 'intcdump.output'.*
* *
************************************************************************
character dumpname*24
************************************************************************
* *
* Logitname is the name to use for the output file which will contain *
* the summary results from the logistic regression fit. This output *
* will contain estimated coefficients with their standard errors and *
* t-values as well as a set of statistics containing the log-likeli- *
* hood of the model, the corresponding value of Bic, total number of *
* woman exposure-years, total number of births and an error return *
* code (which is set to 0 when all has executed normally). The *
* default name is 'logitresults'. *
* *
************************************************************************
character logitname*24
************************************************************************
* *
* Gibbsname is the name to use for the output file which will contain *
* the summary results from the estimated fit resulting from running *
* the Gibbs sampler. Currently this output contains the final *
* estimates of both the fixed and random effect parameters as well as *
* the final estimate of the variance of the random effect parameters *
* and, lastly, the log-likelihood of the full model at these final *
* estimates. The default name is 'gibbsresults'. *
* *
************************************************************************
character gibbsname*24
************************************************************************
* *
* Residname is the name to use for the separate output file created *
* when the residout option is selected. The default name is *
* 'resids.output'. *
* *
************************************************************************
character residname*24
************************************************************************
* *
* Mlfrprioname is the name to use for the output file which will *
* contain the summary results from the marginal likelihood based on *
* draws from the joint prior distribution program, mlfrprio. This *
* output consists of a number of records containing model recognition *
* information, the fixed effects prior distribution covariance hyper- *
* parameter matrix, the marginal likelihood along with its variance *
* and the log marginal likelihood along with its asymptotic variance. *
* The default name is 'mlfrprio.output'. *
* *
************************************************************************
character mlfrprioname*24
************************************************************************
* *
* Printopt is a two character string containing the user's choice of *
* print options. The first character is the number of decimal places *
* the output results should be printed to. The second character is *
* either a 'c' or a ' '. If it is a 'c', the wfsfit program is to *
* print the entire estimated correlation matrix of the parameters as *
* part of its output. If it is a ' ', the matrix is not output. *
* *
************************************************************************
character printopt*2
************************************************************************
* *
* Provset is a character string of up to 62 characters consisting of *
* a comma separated list of expressions of either the form *
* or -, where may be any integer between 0 *
* and 22, inclusive. The provset option permits the user to specify *
* which province delimited subsets of the sample of Iranian women to *
* use for estimating the current model. Use of the provset option *
* causes either the wfsfit or wfsgibbs program to calculate parameter *
* estimates applicable to specific provinces, or ranges of provinces. *
* The expressions designate individual province numbers of *
* provinces to be included in the current model estimation. The *
* - expressions designate ranges of province numbers *
* to be included. The first in each range must be less than *
* or equal to the second in the same range expression. The *
* range expressions are inclusive of the province numbers specified. *
* All provinces included in any of the separate expressions will be *
* used during the current model estimation. Extra blanks before or *
* after the delimiters (',' and '-') will be ignored. The default *
* value is '0-22'. Other sample valid provset options are: '0-1,20', *
* '9-10,12-20', '0-1,5,7-10,1,20-22,19' or '0', and so forth. *
* *
************************************************************************
character provset*62
************************************************************************
* *
* Fixedmethod is a character string identifying which method is to be *
* used for simulating from the univariate posterior distribution of *
* each fixed effect parameter estimate. The default (and currently *
* only) value for the fixedmethod option is 'metropolis'. *
* *
************************************************************************
character fixedmethod*24
************************************************************************
* *
* Fixedprior is a character string identifying the form of the prior *
* distribution of the fixed effect parameter estimates. The range of *
* the support for this prior distribution may be specified using the *
* fixedsupport option. The default value for the fixedprior option *
* is 'normal'. The only other alternative currently implemented is *
* 'uniform'. *
* *
************************************************************************
character fixedprior*24
************************************************************************
* *
* Fixedsupport is a character string specifying the lower and upper *
* bounds of the support of the prior distribution of the fixed effect *
* parameter estimates. Fixedsupport must consist of two integers *
* (or the special value 'inf') separated by a ':'. The first integer *
* should be preceded by a minus sign (but no error will be flagged if *
* this is not the case). An example value for fixed support might be *
* '-1000:1000'. The default support bounds option is '-inf:inf'. *
* *
************************************************************************
character fixedsupport*24
************************************************************************
* *
* Randprior is a character string identifying the form of the prior *
* distribution of the random effect parameter estimates. The default *
* (and currently only) value for the randprior option is 'normal'. *
* *
************************************************************************
character randprior*24
************************************************************************
* *
* Varprior is a character string identifying the form of the prior *
* distribution of the variance of the random effect parameter *
* estimates. The default (and currently only) value for the varprior *
* option is 'inverted gamma'. *
* *
************************************************************************
character varprior*24
************************************************************************
* *
* Fixedmargvar is a double precision number containing the multi- *
* variate normal prior distribution variance for all the fixed effect *
* parameter estimates, with the exception of the intercept, where the *
* covariates have been standardized to have mean 0 and variance 1. *
* It is not necessary for the covariates to actually be standardized. *
* See University of Washington Department of Statistics technical *
* report 255, "Approximate Bayes factors and accounting for model *
* uncertainty in generalized linear models", by Adrian Raftery, for a *
* discussion regarding the use of reference multivariate normal prior *
* distributions as approximations to improper prior distributions in *
* the situation where there is little prior information. The default *
* value for fixedmargvar is e=2.7182818d0. *
* *
************************************************************************
double precision fixedmargvar
************************************************************************
* *
* There are two parameters available in the specification of the *
* distribution of the variance of the random effect parameter *
* estimates. These are a shape parameter, varshape, and a scale *
* parameter, varscale. For the prior distribution to be a proper *
* distribution both of these parameters should be positive numbers. *
* This warning definitely applies using the default form for the *
* prior distribution, namely the 'inverted gamma'. However it is *
* worth observing that by setting both of these parameters to zero, *
* the 'inverted gamma' distribution reduces down to a Jeffreys prior, *
* which if used leads to an improper posterior distribution. The *
* default value for varshape is 0.47d0 and for varscale is 0.0226d0. *
* *
************************************************************************
double precision varshape
double precision varscale
************************************************************************
* *
* Firstvariance is a double precision number containing a number to *
* use as the initial value for the variance of the random effects in *
* the wfsgibbs program. So long as the sampler is mixing well, the *
* precise setting of this option should not be terribly important. *
* The default value for firstvariance is 0.25d0. *
* *
************************************************************************
double precision firstvariance
************************************************************************
* *
* Iteracnt is the total number of complete iterations of the Gibbs *
* sampler to be executed. One complete iteration is a single pass *
* through all of the parameters being estimated. The default number *
* of iterations is 75500. *
* *
************************************************************************
integer iteracnt
************************************************************************
* *
* Burnin is the number of initial complete iterations of the Gibbs *
* sampler which are to be bypassed before calculating any parameter *
* estimates. It is the number of iterations to run the sampler so *
* that the sampled values will (with high probability) be coming from *
* the stationary distribution of the Markov chain. As above, one *
* complete iteration is a single pass through all of the parameters *
* being estimated. The default number of burnin iterations is 500. *
* *
************************************************************************
integer burnin
************************************************************************
* *
* Priosize is the number of draws which the marginal likelihood from *
* the joint prior distribution program, mlfrprio, is to base its *
* calculation on. The default sample size is 1 million. *
* *
************************************************************************
integer priosize
************************************************************************
* *
* Parmdps is the number of places past the decimal point to which the *
* parameter estimates calculated by wfsfit and/or wfsgibbs are to be *
* output. This option is not currently implemented. The default *
* number of decimal places is 6. *
* *
************************************************************************
integer parmdps
************************************************************************
* *
* Cseed and tseed are the 2 seed values required to use the random *
* number functions in the rans library available on Sun Unix systems. *
* For details on how these routines are to be invoked, including any *
* requirements regarding the seed values, see the man pages on rans, *
* setsd, getsd, etc. The default value for cseed is 43827 and the *
* default value for tseed is 1073. *
* *
************************************************************************
integer cseed
integer tseed
************************************************************************
* *
* Randtracecnt is the number of random effect parameter estimates to *
* write to each record of the randomtrace file. Since one record is *
* written per iteration, randtracecnt is the number of random effect *
* parameters which will be available for further analysis after a run *
* of the wfsgibbs program. The random effect parameter estimates *
* traced will be the randtracecnt women with the lowest woman identi- *
* fication numbers among women actually included in the sample being *
* used to estimate the current model. The randtracecnt option must *
* be less than or equal to 1000. The default value for randtracecnt *
* is to trace the first 15 random effect parameters only. *
* *
************************************************************************
integer randtracecnt
************************************************************************
* *
* Fixedcnt is the number of fixed effect parameter realizations from *
* a single univariate posterior distribution to be returned each time *
* the feffect subroutine is called. The default is 1 realization per *
* subroutine call. *
* *
************************************************************************
integer fixedcnt
************************************************************************
* *
* Randcnt is the number of random effect parameter realizations from *
* a single univariate posterior distribution to be returned each time *
* the reffect subroutine is called. The default is 1 realization per *
* subroutine call. *
* *
************************************************************************
integer randcnt
************************************************************************
* *
* Fixedndx is the fixed effect paramter subscript of the fixed effect *
* to be simulated from in the ucfront program. In other words, it is *
* the column subscript into intcovar of the covariate whose univari- *
* ate conditional distribution is to be simulated from. For historic *
* reasons only, the default value is 2. *
* *
************************************************************************
integer fixedndx
************************************************************************
* *
* Womanid is the woman identification number of the random effect to *
* be simulated from in the ucfront program. In other words, only *
* those rows in intcovar associated with this woman are to be used *
* while simulating from the univariate conditional distribution. The *
* default woman id number is 3911 (the woman with the most births). *
* *
************************************************************************
integer womanid
************************************************************************
* *
* Ieee is a single character indicating whether or not the ieee *
* floating point exception handling routines, handlon and handloff, *
* are to be called. 'T' means these routines should be called and *
* 'F' means they should not. 'T' is the default value. *
* *
************************************************************************
character ieee*1
************************************************************************
* *
* Standardize is a single character indicating whether or not a *
* correlation transformation (such as described by Neter, Wasserman *
* and Kutner 1989) is to be applied to the internal covariate matrix. *
* The first column, which must contain an intercept, is never trans- *
* formed. 'T' indicates that the transformation should be applied to *
* the matrix and 'F' means it should not. The default is 'F'. *
* *
************************************************************************
character standardize*1
************************************************************************
* *
* Intcdump is a single character indicating whether or not a dumpout *
* of the internal covariate matrix, intcovar, is to be written to a *
* separate output file. The name to use for this output file may be *
* specified using the dumpname option. 'T' indicates that the dump *
* is to be produced. 'F' indicates that no dump is to be produced. *
* The default is 'F'. *
* *
************************************************************************
character intcdump*1
************************************************************************
* *
* Residout is a single character indicating whether or not a separate *
* output file of fitted value and residuals is to be written out. *
* The name to use for this output file may be specified using the *
* residname option. 'T' indicates that the separate output file is *
* to be written. 'F' indicates that this file should not be written. *
* The default is 'F'. *
* *
************************************************************************
character residout*1
************************************************************************
* *
* Gibbstrace is a single character indicating whether or not the *
* Gibbs sampler trace files (currently using names of 'fixedtrace', *
* 'randomtrace' and 'variancetrace') should be written out by the *
* wfsgibbs program. 'T' indicates the trace files should be output. *
* 'F' indicates they should not be written. The default is 'T'. *
* *
************************************************************************
character gibbstrace*1
************************************************************************
* *
* Gibbstime is a single character indicating whether or not elapsed *
* computer time information used for simulating from the univariate *
* posterior distributions of the parameters should be kept track of *
* by the wfsgibbs program. If so, the total elapsed computer time *
* spent simulating from the posterior distribution of fixed effect *
* parameters, from the posterior distribution of random effect *
* parameters and from the posterior distribution of random effect *
* variance component parameters will be written to separate lines in *
* the wfsgibbs timing summary output file, currently named *
* 'time.summary'. 'T' indicates that this timing information should *
* be output. 'F' indicates that this timing information should not *
* be written. The default is 'T'. *
* *
************************************************************************
character gibbstime*1
************************************************************************
* *
* Mlfrpriotrace is a single character indicating whether or not the *
* marginal likelihood from the joint prior distribution program, *
* mlfrprio, is to include a trace of its calculation of the marginal *
* likelihood in its output file or not. 'T' indicates that the trace *
* should be included in the output. 'F' indicates that the trace *
* should not be included. The default is 'F'. *
* *
************************************************************************
character mlfrpriotrace*1
************************************************************************
* *
* The following 3 selection flags are used by the ucfront check out *
* program. They may be set to specify which univariate conditional *
* distribution(s) are to be simulated from and which are not to be *
* simulated from. 'T' indicates the corresponding distribution *
* (fixed = fixed effect parameters, rand = random effect parameters *
* and var = variance of the random effect parameters) is to be simu- *
* lated from. 'F' indicates the corresponding distribution is not to *
* be simulated from. The default for all three is 'F'. *
* *
************************************************************************
character fixedchk*1
character randchk*1
character varchk*1
************************************************************************
* *
* In order to include any of the province level variables we need to *
* know the province number of the province in which the woman lives. *
* There is only one woman missing her province number and not missing *
* any of the other requisite covariates. This is woman 4707. The *
* keep4707 is a "temporary" option which permits me to run models not *
* involving any province level variables without excluding woman 4707 *
* from the sample. If keep4707 is 'T', the chkevent subroutine will *
* not reject intervals belonging to woman 4707 despite the fact that *
* her province number is missing. The chkevent subroutine does this *
* by assuming that woman 4707 lives in province 16 (since women 4706 *
* and 4708 live in province 16). This option permits estimation of *
* previously calculated models without decreasing the sample size. *
* The default setting for this option is 'F'. *
* *
************************************************************************
character keep4707*1
************************************************************************
* *
* When the Iran Fertility Survey programs were originally written, a *
* conservative upper limit of 24 years was placed on all birth *
* intervals. If the input event history file included any interval *
* longer than this, it was assumed to be in error and the interval *
* was excluded from the estimation of the model. This upper limit is *
* certainly reasonable for any closed interval, but it would not be *
* at all surprising for an open interval to be longer than this. The *
* longeropen option was added as another "temporary" option. If this *
* option is set to 'F', the chkevent subroutine executes as it did in *
* the original implementation. If this option is set to 'T', then *
* open intervals of up to 38 years in length are deemed by chkevent *
* to be acceptable and as such are left in the estimation of the *
* model. The default setting for this option is 'T'. *
* *
************************************************************************
character longeropen*1
SHAR_EOF
chmod 644 'optioncb.h'
fi
exit 0
# End of shell archive