"nevtot: total number of events for all individuals"
CALC nevtot = nval(indiv)
CALC nevtot2 = nevtot - 1
"nindivs: number of individuals"
CALC nindivs = nlevels(indiv)
"lex: number of risk group cut points"
CALC lex = nval(ex)
"lhist: length of history variate (all cut points for 1 individual)"
CALC lhist = lex + nval(t) + 2
"lvars: number of intervals for 1 individual (inc. 0 intervals)"
CALC lvars = lhist - 1
"lall: number of intervals for all individuals"
CALC lall = nindivs*lvars
"z: non-zero if the first event for each individual"
VARI [nvalues = nevtot] z
CALC z$[1] = 1
CALC z$[2...nevtot] = indiv$[2...nevtot] - indiv$[1...nevtot2]
"restrict sta, end, ex to the first record only"
SUBSET [condition = z .NE. 0] oldvector=sta,end,ex[]
"if covariates exist, restrict to first record and set up variates 'cov'"
CALC ncovs = nval(covariates)
IF ncovs .GT. 0
SUBSET [condition = z .NE. 0] oldvector=covariates[]
VARI [nvalues = lall] cov[1...ncovs]
ENDIF
"set up variates for the following:\
individual: factor for each individual\
exposure_group: exposure group risk factor\
age_group: age group factor\
interval: interval length\
nevents: no. of events\
(1st 3 will be converted to factors before model fitting)"
VARI [nvalues = lall] individual,exposure_group,age_group,interval,nevents
DELETE [redefine = yes] nevtot,nevtot2,z,lall
"loop for each individual"
FOR ind = 1...nindivs;
"pos1 & pos2: where to place data for each individual"
CALC pos1 = ind * lvars - lvars + 1
CALC pos2 = ind * lvars
"factor for each individual"
CALC individual$[pos1...pos2] = ind
"covariates - if they exist"
IF ncovs .GT. 0
FOR i = 1...ncovs
CALC cov[i]$[pos1...pos2] = covariates[i]$[ind]
ENDFOR
ENDIF
"exposure risk group cut points"
VARI [nvalues = lex] excuts
CALC excuts$[1...lex] = ex[]$[ind]
"hist: ordered cut points (start, end, age and risk group)"
POINTER [VALUES=(sta, end)$[ind], excuts, t] histp
VARI [nvalues = lhist] hist
EQUATE histp; newstructures = hist
SORT hist
"number of events in each interval"
SUBSET [condition = indiv .EQ. ind] oldvector = eventday;\
newvector = eventdayi
TALLY [print=*] eventdayi; limits = hist+0.5; frequencies = neventsi
CALC nevents$[pos1...pos2] = neventsi$[2...lhist]
"exposure group risk factor"
GROUPS [lmethod = give] vector = hist; factor = exi;\
limits = excuts+0.5; levels = exc
CALC exposure_group$[pos1...pos2] = exi$[2...lhist]
"age group factor"
GROUPS [lmethod = *] vector = hist; factor = agegri; limits = t+0.5
CALC age_group$[pos1...pos2] = agegri$[2...lhist]
"interval lengths"
CALC interval$[pos1...pos2] = hist$[2...lhist] - hist$[1...lvars]
"set interval = 0 if cut points are less than sta or greater\
than end (these will be taken out before fitting the model)."
FOR i=1...lvars; k=pos1...pos2
IF hist$[i] .LT. sta$[ind]
CALC interval$[k] = 0
ENDIF
ENDFOR
FOR j=2...lhist; k=pos1...pos2
IF hist$[j] .GT. end$[ind]
CALC interval$[k] = 0
ENDIF
ENDFOR
ENDFOR
DELETE [redefine = yes] covariates,ex,lex,histp,eventday,eventdayi,indiv,\
neventsi,exi,excuts,exc,agegri,t,hist,sta,end,pos1,pos2,lhist,lvars
"cut out all intervals in which the interval length is 0 and\
convert covariates to factors"
IF ncovs .GT. 0
RESTRICT individual,interval,age_group,exposure_group,nevents,cov[];\
condition = interval .NE. 0
GROUPS [redefine=yes] cov[]
ELSE
RESTRICT individual,interval,age_group,exposure_group,nevents;\
condition = interval .NE. 0
ENDIF
"log the intervals"
CALC loginterval = LOG(interval)
"convert individual, agegr and vacc to factors"
GROUPS [redefine=yes] individual,age_group,exposure_group
"... and finally fit the model :-)"
MODEL [distribution = poisson; link = logarithm;\
offset = loginterval; groups = individual] nevents
FIT [print = model, summary] ##fit
RKEEP Y = nevents; estimates = estimate; se = s_e; lower = CI_low;\
upper = CI_hi
CALC RR = exp(estimate)
CALC CI_low = exp(CI_low)
CALC CI_hi = exp(CI_hi)
"print out estimates, SEs, RRs and CIs for RRs"
PRINT estimate, s_e, RR, CI_low, CI_hi; decimals = 3
"table of total number of events by age group and vaccine risk factor"
TABULATE [print = totals; classification = age_group,exposure_group] nevents