/******************************************************************************
Project: Immigration Data Linkage Project (ircc_link)

File: S:\ircc_link\prog\gillesd\concept_samples\first_hc_contact.sas

Programmers: Gilles R. Detillieux and Heather Prior

First Created: Oct. 20, 2022

Revisions: none

Input datasets:
registry.ircc_crosswalk_1985jan: IRCC linkage dataset
registry.ircc_legacy_landing_1985jan: IRCC landing data from 1985 to 1999
registry.ircc_landing_2000jan: IRCC landing data from 2000 to 2017
registry.mhmrs_1970jun: MH insurance registry
health.mhhosp_*: MH hospital abstracts
health.mhmed_*: MH medical services claims

Description: This program finds landing date of immigrants living in Manitoba,
for anyone listed in the IRCC immigration landing datasets and crosswalk,
then gets the first MB health insurance registry coverage episode ending after
landing. A cohort is formed from immigrants and all other Manitobans in the
study period, to do a comparative analysis of first contact with health care
between the two groups. Health care contact can be either an in-patient
hospitalization of an ambulatory visit.

******************************************************************************/

options nocenter;


*** Get immigrant cohort from IRCC crosswalk file, include all valid linkages;
*** and sort by IRCC's IDNO identifier;
proc sort data=registry.ircc_crosswalk_1985jan
                 (keep=filephin filephintype idno acqdt
                  where=(filephintype in ('0','1') and acqdt = '14may2019'd))
          out=ircccw (drop=filephintype rename=(filephin=scrphin)) nodupkey;
   by idno filephin;
run;

*** Get arrival & landing dates from legacy and current landing datasets;
data landing;
   set registry.ircc_legacy_landing_1985jan
         (in=in_ll keep=idno arrival_dt landing_dt acqdt)
       registry.ircc_landing_2000jan
         (in=in_cl keep=idno arrival_dt landing_dt acqdt);
   where acqdt = '14may2019'd;
run;

*** sort landing data, then merge with crosswalk data, by IDNO identifier;
proc sort data=landing;
   by idno landing_dt;
run;

data landing;
   merge ircccw (in=in_cw) landing (in=in_la);
   by idno;
   * keep all individuals with linkage and landing records;
   if in_cw & in_la;
run;

*** sort by PHIN for duplicate detection;
proc sort data=landing;
   by scrphin landing_dt;
run;

*** avoid duplicate PHINs and keep earliest landing record per person;
data landing;
   set landing;
   by scrphin landing_dt;
   * keep one record per valid PHIN;
   if first.scrphin & scrphin ^= .;
run;


*** coverage episodes pulled from registry using regcov macro, for all MB;
*** ignore coverage gaps of 360 days or less;
%regcov(dsname=,
        startdt=01jan1985,
        enddt=30nov2019,
        id=scrphin,
        outds=regcov,
        gapdays=360);

*** sort by person so we can merge with grouped landing data;
proc sort data=regcov;
   by scrphin covdt enddt;
run;

*** merge in landing data to find coverage episodes during or after landing;
*** keep a cohort of landed immigrants and other Manitoba residents;
*** with coverage starting within our study period;
data cohort;
   merge landing (in=in_la)
         regcov (in=in_co);
   by scrphin;
   imm = in_la;
   if in_co & covdt >= '01jan1985'd & (^imm | enddt >= landing_dt);
   label imm = 'Person is a landed immigrant';
run;

*** keep only first coverage episode on or after landing;
*** or first episode within study period for other Manitobans;
data cohort;
   set cohort;
   by scrphin;
   if first.scrphin;
   *** drop dupl. or registered in error (may show up under another phin);
   if canccode in ('W','5','6') then delete;
   pop = 1;
   startdt = max(covdt, landing_dt);
   followupdt = min(enddt, '31dec2019'd);
   * number of person-years from later of start of coverage or landing date;
   * to earlier of end of coverage or end of study period;
   pyears = yrdif(startdt, followupdt+1, 'actual');
   if pyears < 0 then delete;
   label pyears = 'Number of person-years as resident of Manitoba';
run;


*** make a format from list of scrambled PHINs, for selecting health records;
data control;
   set cohort (keep=scrphin) end=eof;
   retain fmtname 'scrphinf';
   rename scrphin=start;
   label = 'Y';
   output;
   if eof then do;
      HLO = 'O';
      label = 'N';
      scrphin = 0;
      output;
   end;
run;

proc format cntlin=control;
run;


*** get hospital abstracts;
data hosp;
   set health.mhhosp_1983apr (in=haum keep=scrphin acqdt admdt sepdt
                                           hosp transact abstype)
       health.mhhosp_1987apr (in=haum keep=scrphin acqdt admdt sepdt
                                           hosp transact abstype)
       health.mhhosp_1995apr (in=haum keep=scrphin acqdt admdt sepdt
                                           hosp transact abstype)
       health.mhhosp_1997apr (in=haum keep=scrphin acqdt admdt sepdt
                                           hosp transact abstype)
       health.mhhosp_2004apr_rec (in=haum keep=scrphin acqdt admdt sepdt
                                           hosp transact abstype)
       health.mhhosp_2004apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2007apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2008apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2009apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2011apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2012apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2013apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2014apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2015apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2016apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2017apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2018apr_dad (in=dad keep=scrphin acqdt provcd admdt sepdt
                       hosp transact abstracttype rename=(abstracttype=abstype))
       health.mhhosp_2008apr_rec (in=oop10 keep=scrphin acqdt admdt sepdt
                                           hosp transact abstype)
       health.mhhosp_2012apr_rec (in=oop10 keep=scrphin acqdt admdt sepdt
                                           hosp transact abstype)
       health.mhhosp_2009apr_recicd9 (in=oop9 keep=scrphin acqdt admdt sepdt
                                           hosp transact abstype);

    where put(scrphin, scrphinf.) = 'Y' &
         '01jan1985'd <= acqdt <= '09jan2020'd &
         '01jan1985'd <= sepdt <= '31mar2019'd &
         transact = '1' &
         abstype ^= '4';

   if dad & provcd ^= 'MB' then delete;
   drop provcd;

   if dad | oop10 then icd = '10';
   else icd = '9';

   *** exclude non-acute care facilities;
   hosptype = put(hosp, $hsptypf.);
   if hosptype in ('07','08','09','10') then delete;
run;

*** remove duplicates;
proc sort data = hosp;
   by scrphin admdt sepdt hosp descending acqdt;
run;

data hosp;
   set hosp;
   by scrphin admdt sepdt hosp descending acqdt;
   if first.hosp;
run;

proc sort data = hosp nodupkey;
   by scrphin sepdt;
run;


*** get medical services claims;
data medin;
   set health.mhmed_1973mar (keep=scrphin acqdt servdt diag prefix tariff
                                  md refmd mdbloc nserv fee opd ngc)
       health.mhmed_1983apr (keep=scrphin acqdt servdt diag prefix tariff
                                  md refmd mdbloc nserv fee opd ngc)
       health.mhmed_1989apr (keep=scrphin acqdt servdt diag prefix tariff
                                  md refmd mdbloc nserv fee opd ngc)
       health.mhmed_1994apr (keep=scrphin acqdt servdt diag prefix tariff
                                  md refmd mdbloc nserv fee opd ngc)
       health.mhmed_1997apr (keep=scrphin acqdt servdt diag prefix tariff
                                  md refmd mdbloc nserv fee opd ngc mgc)
       health.mhmed_2018sep (keep=scrphin acqdt servdt diag prefix tariff
                        md_code refmd_code mdbloc nserv fee opd ngc mgc);
   where '01jan1985'd <= acqdt <= '09jan2020'd &
         '01jan1985'd <= servdt <= '31mar2019'd;
   if prefix = '7' & put(scrphin, scrphinf.) = 'Y';
run;

*** run the net_med macro (DONT exclude nurse practitioners, MDBLOC=200);
%net_med(meddata=medin,
         hospdata=hosp,
         out=med,
         exclude_mdbloc='053' '201' '113' '99',
         type=ambvis,
         debug=debug);

proc delete data=control medin;
run;

*** keep ambulatory visits, remove duplicates;
data med;
   set med (keep=scrphin servdt ambvis);
   where ambvis = 1;
   drop ambvis;
run;

proc sort data = med nodupkey;
   by scrphin servdt;
run;

*** combine both sets of health care usage records;
data hcu;
   merge hosp(in=in_h keep=scrphin sepdt rename=(sepdt=hcudt))
         med(in=in_m keep=scrphin servdt rename=(servdt=hcudt));
   by scrphin hcudt;
   in_hsp = in_h;
   in_med = in_m;
run;

proc delete data=hosp med;
run;

*** merge cohort & HCU, keep records in each person's coverage period;
data hcu;
   merge cohort (in=in_co keep=scrphin startdt followupdt)
         hcu (in=in_hc);
   by scrphin;
   if in_co & in_hc & startdt <= hcudt <= followupdt;
run;

*** sort by date and keep first health care contact;
proc sort data = hcu;
   by scrphin hcudt;
run;

data hcu;
   set hcu;
   by scrphin hcudt;
   if first.scrphin;
run;

proc freq data = hcu;
   tables in_hsp*in_med /list missing;
run;

*** merge first contact back into cohort data;
data final;
   merge cohort (in=in_co)
         hcu (in=in_hc keep=scrphin hcudt);
   by scrphin;
   if in_co;
   event = in_hc;
   censor = 1 - event;
   years2hcu = yrdif(startdt, min(followupdt, hcudt), 'actual');
   days2hcu = min(followupdt, hcudt) - startdt;
run;

*** run a few test statistics;
proc freq data=final;
   tables censor*event / list missing;
run;

proc means data=final n nmiss mean median min max;
   class imm event;
   var years2hcu days2hcu;
run;

*** run a Cox proportional hazard regression, calculate confidence intervals;
proc phreg data=final;
   ods output ParameterEstimates=hr;
   class imm (ref='0') / param=ref;
   model days2hcu*censor(1) = imm;
   title 'Unadjusted PHReg: Time to First Healthcare Contact';
run;

data hr;
   set hr;
   lcl_est = estimate - 1.96*stderr;
   ucl_est = estimate + 1.96*stderr;
   lcl_hr = exp(lcl_est);
   ucl_hr = exp(ucl_est);
run;


*** final output statistics;
proc means data=final n mean stddev min q1 median q3 max;
   class imm;
   var days2hcu;
   title 'Mean & Median # Days from Landing Date to First Healthcare Contact for Immigrants vs. Inter-provincial Migrants';
   title2 'Years of data included: January 1, 1985 to March 31, 2019';
   footnote2 '(First Healthcare Contact could be Ambulatory Visit or Inpatient Hospitalization)';
run;

proc print data=hr noobs;
   var parameter classval0 hazardratio lcl_hr ucl_hr estimate
       stderr chisq probchisq;
   title 'Unadjusted Hazard Ratios of Time to First Healthcare Contact for Immigrants vs. Inter-provincial Migrants';
   title2 'Years of data included: January 1, 1985 to March 31, 2019';
run;

proc univariate data=final;
   class imm;
   var years2hcu;
   histogram;
   title 'Distribution of Years to First Healthcare Contact for Immigrants vs. Inter-provincial Migrants';
   title2 'Years of data included: January 1, 1985 to March 31, 2019';
run;

proc lifetest data=final method=km notable plots=(s) graphics cs=none;
   time years2hcu*censor(1);
   strata imm;
   title 'Kaplan Meier Curve of Years to First Healthcare Contact for Immigrants vs. Inter-provincial Migrants';
run;


*** a little cleanup;
proc delete data=ircccw landing regcov cohort hcu final hr;
run;