Notebook 1 of 5: Identification of nPOD Donors and Matching Process¶
The entire nPOD biobank collection was screened for eligible study donors and then matched in trios, as described in the main text. Each step and the corresponding code is provided below.
Dependencies¶
Set the file location for created data files.
In [1]:
libname files "/mnt/hgfs/myfolders"
Out[1]:
The following macros were used to perform the 1:1:1 matching. The Pharmacoepidemiology Toolbox was developed by Rassen JA and colleagues, is documented and downloadable here. Running the below code chunk for the first time will create a classpath warning message that can be ignored. Re-running the code will resolve this particular error, but not neccessary for execution of this notebook.
In [2]:
%INCLUDE "macros/java_utils.sas";
%INCLUDE "macros/matching.sas";
%toolbox_Start(%QUOTE (macros/pharmacoepi.jar));
PROC javainfo classpaths;
run;
Out[2]:
Data Import & Coding of Needed Variables¶
Obtain all nPOD case records with with donor timeline data from nPOD DataShare.
In [3]:
PROC import out = icu datafile =
"Donor_TimeLine.xlsx"
dbms = xlsx replace;
options validvarname = V7;
run;
Out[3]:
Flags represent potential issues with the timeline and are therefore exlcuded from consideration.
In [4]:
DATA icu1 (drop = npod_caseid);
set icu;
ptid = npod_caseid * 1;
if Admit_Datetime_Flag = 1 then
Hospitalization_Stay_Minutes = .;
if Xclamp_Datetime_Flag = 1 then
Hospitalization_Stay_Minutes = .;
if Xclamp_Datetime_Flag = 1 then
Transport_Duration_Minutes = .;
if Organ_Received_Datetime_Flag = 1 then
Transport_Duration_Minutes = .;
if Organ_Received_Datetime_Flag = 1 then
Processing_Time_Minutes = .;
if Organ_Frozen_Datetime_Flag = 1 then
Processing_Time_Minutes = .;
Hospitalization_Stay_Days = Hospitalization_Stay_minutes / 1440;
Processing_Time_Minutes1 = Processing_Time_Minutes * 1;
transport_duration_hours=Transport_Duration_Minutes / 60;
run;
Out[4]:
Download demographics data from nPOD DataShare. An ERROR message may appear, but safe to ignore.
In [5]:
PROC import datafile = 'demographics.tsv'
out = demographics
dbms = dlm replace;
delimiter = '09'x;
getnames = yes;
guessingrows = 1000;
run;
quit;
Out[5]:
Join timeline and demographics data tables together
In [6]:
PROC sort data = icu1;
by ptid;
run;
PROC sort data = demographics;
by ptid;
run;
DATA all;
merge icu1 (in = a) demographics (in = b);
by ptid;
if a;
run;
Out[6]:
Explore hospitalization stay time (in days). It can be seen below that values are normally distributed across all nPOD samples.
In [7]:
PROC univariate data = ALL noprint;
histogram Hospitalization_Stay_minutes / kernel
(k = normal
c = 0.8
w = 2.5
color = green)
normal (mu = est
sigma = est
color = red
w = 2.5 )
midpoints = 0 to 40000 by 60;
run;
PROC univariate data = ALL noprint;
histogram Hospitalization_Stay_Days / kernel
(k = normal
c = 0.8
w = 2.5
color = green)
normal (mu = est
sigma = est
color = red
w = 2.5)
midpoints = 0 to 30 by 1;
run;
Out[7]:
The UNIVARIATE Procedure
The UNIVARIATE Procedure
Fitted Normal Distribution for Hospitalization_Stay_Minutes (Hospitalization Stay Minutes)
| Parameters for Normal Distribution | ||
|---|---|---|
| Parameter | Symbol | Estimate |
| Mean | Mu | 6788.519 |
| Std Dev | Sigma | 5509.989 |
| Goodness-of-Fit Tests for Normal Distribution | ||||
|---|---|---|---|---|
| Test | Statistic | p Value | ||
| Kolmogorov-Smirnov | D | 0.1810677 | Pr > D | <0.010 |
| Cramer-von Mises | W-Sq | 2.9069783 | Pr > W-Sq | <0.005 |
| Anderson-Darling | A-Sq | 16.1773878 | Pr > A-Sq | <0.005 |
| Quantiles for Normal Distribution | ||
|---|---|---|
| Percent | Quantile | |
| Observed | Estimated | |
| 1.0 | 686.000 | -6029.633 |
| 5.0 | 1788.000 | -2274.607 |
| 10.0 | 2450.000 | -272.816 |
| 25.0 | 3527.000 | 3072.088 |
| 50.0 | 5206.500 | 6788.519 |
| 75.0 | 8440.000 | 10504.951 |
| 90.0 | 13363.000 | 13849.855 |
| 95.0 | 16117.000 | 15851.646 |
| 99.0 | 31266.000 | 19606.672 |
The UNIVARIATE Procedure
The UNIVARIATE Procedure
Fitted Normal Distribution for Hospitalization_Stay_Days
| Parameters for Normal Distribution | ||
|---|---|---|
| Parameter | Symbol | Estimate |
| Mean | Mu | 4.71425 |
| Std Dev | Sigma | 3.826382 |
| Goodness-of-Fit Tests for Normal Distribution | ||||
|---|---|---|---|---|
| Test | Statistic | p Value | ||
| Kolmogorov-Smirnov | D | 0.1810677 | Pr > D | <0.010 |
| Cramer-von Mises | W-Sq | 2.9069783 | Pr > W-Sq | <0.005 |
| Anderson-Darling | A-Sq | 16.1773878 | Pr > A-Sq | <0.005 |
| Quantiles for Normal Distribution | ||
|---|---|---|
| Percent | Quantile | |
| Observed | Estimated | |
| 1.0 | 0.47639 | -4.18725 |
| 5.0 | 1.24167 | -1.57959 |
| 10.0 | 1.70139 | -0.18946 |
| 25.0 | 2.44931 | 2.13339 |
| 50.0 | 3.61563 | 4.71425 |
| 75.0 | 5.86111 | 7.29510 |
| 90.0 | 9.27986 | 9.61795 |
| 95.0 | 11.19236 | 11.00809 |
| 99.0 | 21.71250 | 13.61574 |
Hospitalization stay time (in days) will be binned based on the definition used by In't Veld and colleagues (2010, Diabetes), and consistent with our previous work. Other coding performed to re-format the data.
In [8]:
DATA all1;
set all;
hospbin =.;
if Hospitalization_Stay_Days^ =. and Hospitalization_Stay_Days < 3 then
hospbin = 3;
if Hospitalization_Stay_Days >=3 and Hospitalization_Stay_Days < 6 then
hospbin = 6;
if Hospitalization_Stay_Days >= 6 then
hospbin = 9;
*if Hospitalization_Stay_Days>=9 then hospbin=12;
if donor_type = "Autoab Pos" then
donor = 3;
if donor_type = "T1D" then
donor = 1;
if donor_type = "T1D Medalist" then
donor = 1;
if donor_type = "T2D" then
donor = 2;
if donor_type = "No diabetes" then
donor = 0;
if donor_type = "Gestational diabetes" then
donor = 99;
if donor_type = "Pregnancy" then
donor = 99;
if donor_type = "Other" then
donor = 99;
if donor_type = "Pending" then
donor = 99;
if donor_type = "Autoab Pos by screening" then
donor = 3;
if donor_type = "Preclinical T2D" then
donor = 99;
if donor_type = "Cystic fibrosis" then
donor = 99;
if donor_type = "Fulminant" then
donor = 99;
if donor_type = "Transplant" then
donor = 99;
if donor_type = "Virus infected cell line" then
donor = 99;
if donor_type = "MODY" then
donor = 99;
if donor_type = "CMA" then
donor = 99;
if donor_type = "T2D+Incretin" then
donor = 2;
if donor_type = "TMA" then
donor = 99;
if donor_type = "Gastric Bypass" then
donor = 99;
* coding for sex and ethnicity so as to be used for matching;
sex = .;
if oppc_gender = "Male" then
sex = 1;
if oppc_gender = "Female" then
sex = 2;
ethnicity = .;
if oppc_ethnicity = "Caucasian" then
ethnicity = 1;
if oppc_ethnicity = "Hispanic/Latino" then
ethnicity = 2;
if oppc_ethnicity = "African Am" then
ethnicity = 3;
if oppc_ethnicity = "Asian" then
ethnicity = 4;
if oppc_ethnicity = "American Indian/Alaska Native" then
ethnicity = 5;
run;
Out[8]:
Case Exclusion¶
Limiting donors to only nPOD cases, i.e. excluding donors from other collections stored at nPOD. Calculated age for donor with pre-term births.
In [9]:
DATA all1;
set all1;
if ptid < 6000 then
delete;
if ptid > 6999 then
delete;
*for those donors born pre-term, age was based on numbe of lived days after birth (i.e. number of days after gestational age);
if ptid = 6120 then
oppc_ageR = 1 / 365;
if ptid = 6164 then
oppc_ageR = 13.75 / 365;
if ptid = 6192 then
oppc_ageR = 1 / 365;
if ptid = 6200 then
oppc_ageR = 2.4 / 365;
if ptid = 6201 then
oppc_ageR = 1 / 365;
if ptid = 6210 then
oppc_ageR = 5.48 / 365;
if ptid = 6214 then
oppc_ageR = 5.95 / 365;
run;
Out[9]:
Further limiting donors to only control cases.
In [10]:
DATA all2;
set all1;
* only looking at controls;
if donor^ = 0 then
delete;
if hospbin = . then
delete;
if (Organ_Received_Datetime_Flag^ = 0 and
Organ_Frozen_Datetime_Flag ^= 0 and
Admit_Datetime_Flag^ = 0 and
Braindeath_Datetime_Flag^ = 0 and
clamp_Datetime_Flag^ = 0) then
delete;
run;
Out[10]:
Preparation for Donor Matching¶
For all donors, estimate 3-group propensity score using polytomous logistic regression
In [11]:
PROC logistic data = all2
order = data DESCENDING;
class hospbin sex ethnicity;
model hospbin(ref="3") = oppc_ager sex oppc_bmiR ethnicity /LINK=glogit;
output out = input_3way_cohort_long(KEEP = ptid hospbin prob) PRED = prob;
run;
quit
Out[11]:
If propensity score not generated for a donor (because of unavailable/missing data), then delete case. Calculate logit of propensity score.
In [12]:
DATA input_3way_cohort_long;
set input_3way_cohort_long;
if prob = . then
delete;
logitpr = log(prob / (1 - prob));
drop prob ;
rename logitpr = prob;
run;
Out[12]:
SD of logit is 1.4180436, .2 of that is .2836087, but that caliper width is recommended for 1:1 matching. No studies have been done on optimal caliper width for 1:1:1 matching. We therefore conservatively doubled the 1:1 caliper width allowing .5672174; authors of the Pharmacoepidemiology Toolbox recommend 3 times the 1:1.
In [13]:
PROC means data = input_3way_cohort_long
n
mean
std;
var prob;
run;
Out[13]:
Need to transform the data from one row per donor per level to one row per donor with one variable per level.
In [14]:
PROC transpose data = input_3way_cohort_long OUT = input_3way_cohort
prefix = prob;
by ptid hospbin;
run;
data input_3way_cohort (drop=_name_);
set input_3way_cohort;
ptid1 = put(ptid, 7.);
drop ptid ;
rename ptid1 = ptid;
run;
Out[14]:
A priori to running the matching macros, case 6210 was identified as having some sample quality issues for this experiment, therefore is being removed as a candidate.
In [15]:
DATA input_3way_cohort1;
set input_3way_cohort;
if ptid = 6210 then
delete;
run;
Out[15]:
Macro Matching¶
Use of the matching macro needs the location of the path to a temp directory. This temp directory changes file path from one session to another. The %LET statement stores the session path, and it is then used below when running the macro. If problems arise, an alternative solution is to run %put %sysfunc(getoption(work)); then copy the file path and paste it after the %QUOTE in the macro below (in particular, you need all the numbers and letters that appear after SAS_work).
In [19]:
* Run 1:1:1 nearest neighbor match;
%match_NearestNeighborMatch (input_3way_cohort1,
output_3way_cohort,
%QUOTE(&file_path),
var_patient_id = ptid,
var_exposure = hospbin,
var_ps = prob1 prob2 prob3,
exp_groups = 3 6 9,
caliper = .5672174,
ratio = 1);
DATA match_import1 (keep = pat_id ptid SET_NUM);
set match_import;
ptid = pat_id * 1;
run;
DATA match_import1 (keep = ptid set_num);
set match_import1;
run;
PROC sort data = match_import1;
by ptid;
run;
PROC sort data = all2;
by ptid;
run;
DATA study;
merge match_import1 (in = a) all2 (in = b);
by ptid;
if a;
run;
PROC freq data = study;
tables hospbin;
run;
Out[19]:
Verification of Results¶
Need to determine if there are any differences in age, sex, BMI, and ethnicity between the three groups, as these variables were used for matching. Here, distribution of variables is examined.
In [20]:
PROC univariate data = study;
var oppc_ager
oppc_bmir
Transport_Duration_Minutes
Processing_Time_Minutes1;
run;
PROC univariate data = study noprint;
histogram Transport_Duration_Minutes / kernel
(k = normal
c = 0.8
w = 2.5
color = green)
normal (mu = est
sigma = est
color = red
w = 2.5 )
midpoints = 0 to 2880 by 60;
run;
PROC univariate data = study noprint;
histogram Processing_Time_Minutes1 / kernel
(k = normal
c = 0.8
w = 2.5
color = green)
normal (mu = est
sigma = est
color = red
w = 2.5)
midpoints = 0 to 400 by 5;
run;
Out[20]:
The UNIVARIATE Procedure
Variable: oppc_ageR
| Moments | |||
|---|---|---|---|
| N | 45 | Sum Weights | 45 |
| Mean | 24.1774795 | Sum Observations | 1087.98658 |
| Std Deviation | 16.3125239 | Variance | 266.098436 |
| Skewness | 0.5133685 | Kurtosis | 0.23193307 |
| Uncorrected SS | 38013.1042 | Corrected SS | 11708.3312 |
| Coeff Variation | 67.4699111 | Std Error Mean | 2.43172749 |
| Basic Statistical Measures | |||
|---|---|---|---|
| Location | Variability | ||
| Mean | 24.17748 | Std Deviation | 16.31252 |
| Median | 22.70000 | Variance | 266.09844 |
| Mode | 0.33000 | Range | 71.99342 |
| Interquartile Range | 24.00000 | ||
Note: The mode displayed is the smallest of 6 modes with a count of 2.
| Tests for Location: Mu0=0 | ||||
|---|---|---|---|---|
| Test | Statistic | p Value | ||
| Student's t | t | 9.942512 | Pr > |t| | <.0001 |
| Sign | M | 22.5 | Pr >= |M| | <.0001 |
| Signed Rank | S | 517.5 | Pr >= |S| | <.0001 |
| Quantiles (Definition 5) | |
|---|---|
| Level | Quantile |
| 100% Max | 72.00000000 |
| 99% | 72.00000000 |
| 95% | 50.00000000 |
| 90% | 45.80000000 |
| 75% Q3 | 38.00000000 |
| 50% Median | 22.70000000 |
| 25% Q1 | 14.00000000 |
| 10% | 0.50000000 |
| 5% | 0.33000000 |
| 1% | 0.00657534 |
| 0% Min | 0.00657534 |
| Extreme Observations | |||
|---|---|---|---|
| Lowest | Highest | ||
| Value | Obs | Value | Obs |
| 0.00657534 | 35 | 45.8 | 28 |
| 0.33000000 | 34 | 46.0 | 6 |
| 0.33000000 | 20 | 50.0 | 5 |
| 0.42000000 | 19 | 51.0 | 29 |
| 0.50000000 | 36 | 72.0 | 9 |
The UNIVARIATE Procedure
Variable: oppc_bmiR
| Moments | |||
|---|---|---|---|
| N | 45 | Sum Weights | 45 |
| Mean | 24.06 | Sum Observations | 1082.7 |
| Std Deviation | 5.45316422 | Variance | 29.737 |
| Skewness | 0.26293153 | Kurtosis | -0.796331 |
| Uncorrected SS | 27358.19 | Corrected SS | 1308.428 |
| Coeff Variation | 22.6648554 | Std Error Mean | 0.81290973 |
| Basic Statistical Measures | |||
|---|---|---|---|
| Location | Variability | ||
| Mean | 24.06000 | Std Deviation | 5.45316 |
| Median | 24.20000 | Variance | 29.73700 |
| Mode | 17.10000 | Range | 21.50000 |
| Interquartile Range | 7.90000 | ||
Note: The mode displayed is the smallest of 2 modes with a count of 2.
| Tests for Location: Mu0=0 | ||||
|---|---|---|---|---|
| Test | Statistic | p Value | ||
| Student's t | t | 29.59738 | Pr > |t| | <.0001 |
| Sign | M | 22.5 | Pr >= |M| | <.0001 |
| Signed Rank | S | 517.5 | Pr >= |S| | <.0001 |
| Quantiles (Definition 5) | |
|---|---|
| Level | Quantile |
| 100% Max | 36.0 |
| 99% | 36.0 |
| 95% | 32.7 |
| 90% | 32.2 |
| 75% Q3 | 27.9 |
| 50% Median | 24.2 |
| 25% Q1 | 20.0 |
| 10% | 17.1 |
| 5% | 16.3 |
| 1% | 14.5 |
| 0% Min | 14.5 |
| Extreme Observations | |||
|---|---|---|---|
| Lowest | Highest | ||
| Value | Obs | Value | Obs |
| 14.5 | 35 | 32.2 | 7 |
| 15.7 | 3 | 32.4 | 30 |
| 16.3 | 26 | 32.7 | 13 |
| 17.0 | 36 | 34.0 | 45 |
| 17.1 | 34 | 36.0 | 14 |
The UNIVARIATE Procedure
Variable: Transport_Duration_Minutes (Transport Duration Minutes)
| Moments | |||
|---|---|---|---|
| N | 45 | Sum Weights | 45 |
| Mean | 808.377778 | Sum Observations | 36377 |
| Std Deviation | 493.579195 | Variance | 243620.422 |
| Skewness | 1.09542492 | Kurtosis | 3.65355162 |
| Uncorrected SS | 40125657 | Corrected SS | 10719298.6 |
| Coeff Variation | 61.0579866 | Std Error Mean | 73.5784422 |
| Basic Statistical Measures | |||
|---|---|---|---|
| Location | Variability | ||
| Mean | 808.3778 | Std Deviation | 493.57920 |
| Median | 835.0000 | Variance | 243620 |
| Mode | 835.0000 | Range | 2656 |
| Interquartile Range | 728.00000 | ||
Note: The mode displayed is the smallest of 2 modes with a count of 2.
| Tests for Location: Mu0=0 | ||||
|---|---|---|---|---|
| Test | Statistic | p Value | ||
| Student's t | t | 10.98661 | Pr > |t| | <.0001 |
| Sign | M | 22.5 | Pr >= |M| | <.0001 |
| Signed Rank | S | 517.5 | Pr >= |S| | <.0001 |
| Quantiles (Definition 5) | |
|---|---|
| Level | Quantile |
| 100% Max | 2715 |
| 99% | 2715 |
| 95% | 1464 |
| 90% | 1248 |
| 75% Q3 | 1124 |
| 50% Median | 835 |
| 25% Q1 | 396 |
| 10% | 178 |
| 5% | 105 |
| 1% | 59 |
| 0% Min | 59 |
| Extreme Observations | |||
|---|---|---|---|
| Lowest | Highest | ||
| Value | Obs | Value | Obs |
| 59 | 6 | 1248 | 43 |
| 84 | 2 | 1456 | 35 |
| 105 | 1 | 1464 | 33 |
| 112 | 5 | 1496 | 3 |
| 178 | 9 | 2715 | 45 |
The UNIVARIATE Procedure
Variable: Processing_Time_Minutes1
| Moments | |||
|---|---|---|---|
| N | 42 | Sum Weights | 42 |
| Mean | 94.2619048 | Sum Observations | 3959 |
| Std Deviation | 68.911417 | Variance | 4748.78339 |
| Skewness | 1.55090872 | Kurtosis | 3.03309837 |
| Uncorrected SS | 567883 | Corrected SS | 194700.119 |
| Coeff Variation | 73.1063277 | Std Error Mean | 10.6332625 |
| Basic Statistical Measures | |||
|---|---|---|---|
| Location | Variability | ||
| Mean | 94.26190 | Std Deviation | 68.91142 |
| Median | 77.50000 | Variance | 4749 |
| Mode | 30.00000 | Range | 325.00000 |
| Interquartile Range | 78.00000 | ||
| Tests for Location: Mu0=0 | ||||
|---|---|---|---|---|
| Test | Statistic | p Value | ||
| Student's t | t | 8.864815 | Pr > |t| | <.0001 |
| Sign | M | 21 | Pr >= |M| | <.0001 |
| Signed Rank | S | 451.5 | Pr >= |S| | <.0001 |
| Quantiles (Definition 5) | |
|---|---|
| Level | Quantile |
| 100% Max | 340.0 |
| 99% | 340.0 |
| 95% | 215.0 |
| 90% | 195.0 |
| 75% Q3 | 120.0 |
| 50% Median | 77.5 |
| 25% Q1 | 42.0 |
| 10% | 30.0 |
| 5% | 25.0 |
| 1% | 15.0 |
| 0% Min | 15.0 |
| Extreme Observations | |||
|---|---|---|---|
| Lowest | Highest | ||
| Value | Obs | Value | Obs |
| 15 | 9 | 195 | 25 |
| 20 | 4 | 200 | 29 |
| 25 | 41 | 215 | 13 |
| 30 | 35 | 255 | 12 |
| 30 | 10 | 340 | 21 |
| Missing Values | |||
|---|---|---|---|
| Missing Value |
Count | Percent Of | |
| All Obs | Missing Obs | ||
| . | 3 | 6.67 | 100.00 |
The UNIVARIATE Procedure
The UNIVARIATE Procedure
Fitted Normal Distribution for Transport_Duration_Minutes (Transport Duration Minutes)
| Parameters for Normal Distribution | ||
|---|---|---|
| Parameter | Symbol | Estimate |
| Mean | Mu | 808.3778 |
| Std Dev | Sigma | 493.5792 |
| Goodness-of-Fit Tests for Normal Distribution | ||||
|---|---|---|---|---|
| Test | Statistic | p Value | ||
| Kolmogorov-Smirnov | D | 0.10088584 | Pr > D | >0.150 |
| Cramer-von Mises | W-Sq | 0.07117205 | Pr > W-Sq | >0.250 |
| Anderson-Darling | A-Sq | 0.64243467 | Pr > A-Sq | 0.091 |
| Quantiles for Normal Distribution | ||
|---|---|---|
| Percent | Quantile | |
| Observed | Estimated | |
| 1.0 | 59.0000 | -339.85913 |
| 5.0 | 105.0000 | -3.48775 |
| 10.0 | 178.0000 | 175.83059 |
| 25.0 | 396.0000 | 475.46367 |
| 50.0 | 835.0000 | 808.37778 |
| 75.0 | 1124.0000 | 1141.29189 |
| 90.0 | 1248.0000 | 1440.92497 |
| 95.0 | 1464.0000 | 1620.24331 |
| 99.0 | 2715.0000 | 1956.61469 |
The UNIVARIATE Procedure
The UNIVARIATE Procedure
Fitted Normal Distribution for Processing_Time_Minutes1
| Parameters for Normal Distribution | ||
|---|---|---|
| Parameter | Symbol | Estimate |
| Mean | Mu | 94.2619 |
| Std Dev | Sigma | 68.91142 |
| Goodness-of-Fit Tests for Normal Distribution | ||||
|---|---|---|---|---|
| Test | Statistic | p Value | ||
| Kolmogorov-Smirnov | D | 0.12503077 | Pr > D | 0.096 |
| Cramer-von Mises | W-Sq | 0.22008976 | Pr > W-Sq | <0.005 |
| Anderson-Darling | A-Sq | 1.48109231 | Pr > A-Sq | <0.005 |
| Quantiles for Normal Distribution | ||
|---|---|---|
| Percent | Quantile | |
| Observed | Estimated | |
| 1.0 | 15.0000 | -66.05002 |
| 5.0 | 25.0000 | -19.08729 |
| 10.0 | 30.0000 | 5.94837 |
| 25.0 | 42.0000 | 47.78186 |
| 50.0 | 77.5000 | 94.26190 |
| 75.0 | 120.0000 | 140.74195 |
| 90.0 | 195.0000 | 182.57544 |
| 95.0 | 215.0000 | 207.61110 |
| 99.0 | 340.0000 | 254.57383 |
After matching completed, list of donors was reviewed. It was determined that tissue was not available for 6001 and 6021, so all blocks that contain these samples were removed from consideration of this experiment and matching verification proceeded.
In [21]:
DATA study1;
set study;
if set_num = 2 then
delete;
if set_num = 12 then
delete;
run;
DATA grouping_list (keep = case match_group oppc_ageR);
set study1;
rename ptid=case;
rename set_num=match_group;
run;
Out[21]:
continuous parametric: ANOVA for multi-group comparisons
In [22]:
%macro anova (X);
PROC anova data = study1;
class hospbin;
model &X = hospbin;
run;
%mend anova;
%anova (oppc_ageR);
run;
%anova (oppc_bmiR);
run;
*transport duration was not formally matched, but checked anyway;
%anova (Transport_Duration_Minutes);
run;
Out[22]:
The ANOVA Procedure
| Class Level Information | ||
|---|---|---|
| Class | Levels | Values |
| hospbin | 3 | 3 6 9 |
| Number of Observations Read | 39 |
|---|---|
| Number of Observations Used | 39 |
The ANOVA Procedure
Dependent Variable: oppc_ageR
| Source | DF | Sum of Squares | Mean Square | F Value | Pr > F |
|---|---|---|---|---|---|
| Model | 2 | 90.855584 | 45.427792 | 0.20 | 0.8180 |
| Error | 36 | 8094.286253 | 224.841285 | ||
| Corrected Total | 38 | 8185.141836 |
| R-Square | Coeff Var | Root MSE | oppc_ageR Mean |
|---|---|---|---|
| 0.011100 | 66.91105 | 14.99471 | 22.40991 |
| Source | DF | Anova SS | Mean Square | F Value | Pr > F |
|---|---|---|---|---|---|
| hospbin | 2 | 90.85558353 | 45.42779177 | 0.20 | 0.8180 |
The ANOVA Procedure
| Class Level Information | ||
|---|---|---|
| Class | Levels | Values |
| hospbin | 3 | 3 6 9 |
| Number of Observations Read | 39 |
|---|---|
| Number of Observations Used | 39 |
The ANOVA Procedure
Dependent Variable: oppc_bmiR
| Source | DF | Sum of Squares | Mean Square | F Value | Pr > F |
|---|---|---|---|---|---|
| Model | 2 | 1.626667 | 0.813333 | 0.02 | 0.9760 |
| Error | 36 | 1202.049231 | 33.390256 | ||
| Corrected Total | 38 | 1203.675897 |
| R-Square | Coeff Var | Root MSE | oppc_bmiR Mean |
|---|---|---|---|
| 0.001351 | 24.02033 | 5.778430 | 24.05641 |
| Source | DF | Anova SS | Mean Square | F Value | Pr > F |
|---|---|---|---|---|---|
| hospbin | 2 | 1.62666667 | 0.81333333 | 0.02 | 0.9760 |
The ANOVA Procedure
| Class Level Information | ||
|---|---|---|
| Class | Levels | Values |
| hospbin | 3 | 3 6 9 |
| Number of Observations Read | 39 |
|---|---|
| Number of Observations Used | 39 |
The ANOVA Procedure
Dependent Variable: Transport_Duration_Minutes Transport Duration Minutes
| Source | DF | Sum of Squares | Mean Square | F Value | Pr > F |
|---|---|---|---|---|---|
| Model | 2 | 688024.667 | 344012.333 | 1.45 | 0.2491 |
| Error | 36 | 8569701.231 | 238047.256 | ||
| Corrected Total | 38 | 9257725.897 |
| R-Square | Coeff Var | Root MSE | Transport_Duration_Minutes Mean |
|---|---|---|---|
| 0.074319 | 57.99492 | 487.9009 | 841.2821 |
| Source | DF | Anova SS | Mean Square | F Value | Pr > F |
|---|---|---|---|---|---|
| hospbin | 2 | 688024.6667 | 344012.3333 | 1.45 | 0.2491 |
continuous non-parametric: wilcoxon rank sum
In [23]:
%macro npar1way (X);
PROC npar1way data = study1 wilcoxon;
class hospbin;
var &X;
title "wilcoxon rank sum for &X on group1";
%mend npar1way;
*processing time was not formally matched, but checked anyway;
%npar1way (Processing_Time_Minutes1);
run;
Out[23]:
The NPAR1WAY Procedure
| Wilcoxon Scores (Rank Sums) for Variable Processing_Time_Minutes1 Classified by Variable hospbin |
|||||
|---|---|---|---|---|---|
| hospbin | N | Sum of Scores |
Expected Under H0 |
Std Dev Under H0 |
Mean Score |
| Average scores were used for ties. | |||||
| 3 | 12 | 191.0 | 228.0 | 30.754403 | 15.916667 |
| 9 | 13 | 254.0 | 247.0 | 31.363461 | 19.538462 |
| 6 | 12 | 258.0 | 228.0 | 30.754403 | 21.500000 |
| Kruskal-Wallis Test | |
|---|---|
| Chi-Square | 1.6532 |
| DF | 2 |
| Pr > Chi-Square | 0.4375 |
p-values were calculated for both 2- and multi-class categorical nominal data. Fisher's exact test was used when values from 1 or more cells in a contingency table (2 x 2 or R x C) were less than 5. Pearson's chi square was used where values in all cells exceeded 5. If values from 2 or more cells in any table were equal to 0, p-values were not generated.
In [24]:
%macro freq (X);
PROC freq data = study1;
tables &X * hospbin / measures chisq;
%mend freq;
run;
%macro freq1 (X);
PROC freq data = study1;
tables &X * hospbin / measures fisher;
%mend freq1;
%freq1(ethnicity); run;
%freq1 (sex);
run;
Out[24]:
Additional Grouping¶
An additional factor is created to specifically test whether differences in the analysis are seen based on age. Groupings based on papers by Gregg et al, J Clin Endocrinol Metab, 2012 (examines Ki67/insulin staining by age groups) and Cole TJ, 2009, Adv Exp Med Biol, 639:1-13 (examines organ development by age).
In [25]:
DATA study1;
set study1;
if oppc_ageR<1 then age_group=1; /*developmental and neonates*/
if oppc_ageR>=1 and oppc_ageR<12 then age_group=1; /*child*/
if oppc_ageR>=12 and oppc_ageR<20 then age_group=2; /*adolescent*/
if oppc_ageR>=20 then age_group=3; /*adult*/
run;
PROC sort data=study1;
by hospbin;
run;
PROC freq data=study1;
tables age_group;
run;
PROC freq data=study1;
tables age_group;
by hospbin;
run;
proc print data=study1;
run;
Out[25]:
Need to test for collinearity of age group with hospitalization time group, since both will be included in the model.
In [26]:
%macro freq (X, Y);
PROC freq data=study1;
tables &X*&Y/measures chisq exact;
%mend freq;
/*ORDINAL VS ORDINAL or 2-class nominal vs. ordinal Tau-b if equal r x c, and tau-c if unequal r x c*/
%freq (hospbin,age_group); run;
Out[26]:
Matching Schema Save¶
The final verified matching scheme is saved for later reference.
In [27]:
libname files "/mnt/hgfs/myfolders";
DATA files.grouping_list (keep = case match_group age_group oppc_ageR);
set study1;
rename ptid=case;
rename set_num=match_group;
run;
Out[27]:
Table 1 - Demographics¶
Additional data files needed¶
The study1 data file does not have info on c-peptide. Downloaded that dataset from nPOD DataShare and merged it.
In [28]:
PROC import datafile = 'Diabetes Information.tsv'
out = diabetesinfo
dbms = dlm replace;
delimiter = '09'x;
getnames = yes;
guessingrows = 1000;
run;
quit;
data diabetesinfo1;
set diabetesinfo;
npodcaseID1=npodcaseid*1;
rename nPODCaseID1=ptid;
cpeptide=oppc_c_peptideR*1;
run;
PROC sort data=diabetesinfo1;
by ptid;
run;
PROC sort data=study1;
by ptid;
run;
data study2;
merge study1 (in=a) diabetesinfo1 (in=b);
by ptid;
if a;
run;
proc print data=study2;
run;
Out[28]:
Data on medications adminstered within 24 hours prior to organ recovery. Downloaded from nPOD DataShare and merged it.
In [37]:
PROC import datafile = 'meds.tsv'
out = meds
dbms = dlm replace;
delimiter = '09'x;
getnames = yes;
guessingrows = 1000;
run;
quit;
data meds1 (drop=npodcaseid);
set meds;
npodcaseID1=npodcaseid*1;
rename nPODCaseID1=ptid;
run;
PROC sort data=meds1;
by ptid;
run;
PROC sort data=study2;
by ptid;
run;
data study3;
merge study2 (in=a) meds1 (in=b);
by ptid;
if a;
run;
proc print data=study3;
run;
Out[37]:
Examined med list (other_meds1, 2, and 3) and it needs to be cleaned up so that each medication is placed into a separate field for further analysis. Given the number of records, the dataset was exported and pre-processed manually, then re-imported into the datafile.
In rare cases, dosages were given. This information was not used in the analysis and thus removed as part of the pre-processing. Generic or brand name of medication was used, but not both. Abbreviations were expanded to full medication names, where distinguishable. Blanks were recoded to 'none'. Misspellings were accounted for in the mappings.
In [65]:
PROC export data=study3
outfile="study3.csv"
dbms=csv
replace;
run;
PROC import datafile = 'study4.csv'
out = study4
dbms = csv replace;
getnames = yes;
guessingrows = 1000;
run;
quit;
PROC print data=study4;
run;
Out[65]:
Processing of medications (entered in open text fields) to categories proceeded in a series of steps outlined below.
In [73]:
/********step 1: Imported groupings of previously identified medications into categories
(see Kaddis et al, AJT, 2010; Campbell-Thompson et al, JAMA, 2012) for use in comparison.*/
PROC import datafile='med_list_categories.txt' out=categories dbms=dlm replace;
delimiter='09'x;
getnames=yes;
guessingrows=1000;
run;
quit;
/********step 2: merge categories with med list*/
DATA opentextmerge;
merge study4 categories;
run;
/********step 3: clean up medication list, i.e. change letters to all lowercase, remove white space between words*/
DATA opentextmerge1;
set opentextmerge;
dmed_other1 = lowcase ( dmed_other1 ) ;
dmed_other2 = lowcase ( dmed_other2 ) ;
dmed_other3 = lowcase ( dmed_other3 ) ;
dmed_other4 = lowcase ( dmed_other4 ) ;
dmed_other5 = lowcase ( dmed_other5 ) ;
dmed_other6 = lowcase ( dmed_other6 ) ;
dmed_other7 = lowcase ( dmed_other7 ) ;
dmed_other8 = lowcase ( dmed_other8 ) ;
dmed_other9 = lowcase ( dmed_other9 ) ;
anestetics = lowcase ( anestetics ) ;
antibiotics = lowcase ( antibiotics ) ;
cardiovascular = lowcase ( cardiovascular ) ;
hormonal = lowcase ( hormonal ) ;
others = lowcase ( others ) ;
fluids_electrolytes = lowcase ( fluids_electrolytes ) ;
dmed_other1 = kcompress ( dmed_other1 ) ;
dmed_other2 = kcompress ( dmed_other2 ) ;
dmed_other3 = kcompress ( dmed_other3 ) ;
dmed_other4 = kcompress ( dmed_other4 ) ;
dmed_other5 = kcompress ( dmed_other5 ) ;
dmed_other6 = kcompress ( dmed_other6 ) ;
dmed_other7 = kcompress ( dmed_other7 ) ;
dmed_other8 = kcompress ( dmed_other8 ) ;
dmed_other9 = kcompress ( dmed_other9 ) ;
anestetics = kcompress ( anestetics ) ;
antibiotics = kcompress ( antibiotics ) ;
cardiovascular = kcompress ( cardiovascular ) ;
hormonal = kcompress ( hormonal ) ;
others = kcompress ( others ) ;
fluids_electrolytes = kcompress ( fluids_electrolytes ) ;
run;
/********step 4: double-check categories dataset to make sure that a drug does not appear in more than 1 category*/
PROC iml;
use opentextmerge1;
/*creation of matrix from variables (n by 1 column matrix for each variable*/
read all var {'anestetics'} into anestetics;
read all var {'antibiotics'} into antibiotics;
read all var {'cardiovascular'} into cardiovascular;
read all var {'hormonal'} into hormonal;
read all var {'others'} into others;
read all var {'fluids_electrolytes'} into fluids_electrolytes;
/*matrix initiation*/
b=138; /*number of medications in longest list of (Anestetics, Antibiotics, Cardiovascular, Hormonal, Others, Fluids_Electrolytes)*/
c=6; /*number of categories in list*/
d=138;
/*create macro and DO loop for matrix calculations*/
%macro evaluation;
dmedrecode=j(b,c,0);
do n=1 to b; /*number of records being examined for each variable*/
do p=1 to d; /*max number of different entries in variable list that need to be checked*/
/*coding notes: anestetics=1, antibiotics=2, cardiovascular=3, hormonal=4, others=5, fluids_electrolytes=6*/
if (anestetics[n]=antibiotics[p]) then dmedrecode[n,1]=1;
if (anestetics[n]=cardiovascular[p]) then dmedrecode[n,1]=1;
if (anestetics[n]=hormonal[p]) then dmedrecode[n,1]=1;
if (anestetics[n]=others[p]) then dmedrecode[n,1]=1;
if (anestetics[n]=fluids_electrolytes[p]) then dmedrecode[n,1]=1;
if (antibiotics[n]=anestetics[p]) then dmedrecode[n,2]=1;
if (antibiotics[n]=cardiovascular[p]) then dmedrecode[n,2]=1;
if (antibiotics[n]=hormonal[p]) then dmedrecode[n,2]=1;
if (antibiotics[n]=others[p]) then dmedrecode[n,2]=1;
if (antibiotics[n]=fluids_electrolytes[p]) then dmedrecode[n,2]=1;
if (cardiovascular[n]=anestetics[p]) then dmedrecode[n,3]=1;
if (cardiovascular[n]=antibiotics[p]) then dmedrecode[n,3]=1;
if (cardiovascular[n]=hormonal[p]) then dmedrecode[n,3]=1;
if (cardiovascular[n]=others[p]) then dmedrecode[n,3]=1;
if (cardiovascular[n]=fluids_electrolytes[p]) then dmedrecode[n,3]=1;
if (hormonal[n]=anestetics[p]) then dmedrecode[n,4]=1;
if (hormonal[n]=antibiotics[p]) then dmedrecode[n,4]=1;
if (hormonal[n]=cardiovascular[p]) then dmedrecode[n,4]=1;
if (hormonal[n]=others[p]) then dmedrecode[n,4]=1;
if (hormonal[n]=fluids_electrolytes[p]) then dmedrecode[n,4]=1;
if (others[n]=anestetics[p]) then dmedrecode[n,5]=1;
if (others[n]=antibiotics[p]) then dmedrecode[n,5]=1;
if (others[n]=hormonal[p]) then dmedrecode[n,5]=1;
if (others[n]=cardiovascular[p]) then dmedrecode[n,5]=1;
if (others[n]=fluids_electrolytes[p]) then dmedrecode[n,5]=1;
if (fluids_electrolytes[n]=anestetics[p]) then dmedrecode[n,6]=1;
if (fluids_electrolytes[n]=antibiotics[p]) then dmedrecode[n,6]=1;
if (fluids_electrolytes[n]=hormonal[p]) then dmedrecode[n,6]=1;
if (fluids_electrolytes[n]=cardiovascular[p]) then dmedrecode[n,6]=1;
if (fluids_electrolytes[n]=others[p]) then dmedrecode[n,6]=1;
end;
end;
/*create matrix*/
create compare from dmedrecode;
append from dmedrecode;
close compare;
%mend evaluation;
%evaluation;
quit;
data merged;
merge categories compare;
run;
proc freq data=compare;
tables col1 col2 col3 col4 col5 col6;
run;
/*NOTE: To interpret the results, compare the number of medications in each category to the frequency of 0's in the column.
If they match, medications are not repeated across categories. For example, there are 66 medications listed under anestetics and
the freq of 0's in col1 is 66; thus successful. All columns pass*/
/********step 5: code open text medication fields into categories*/
proc iml;
use opentextmerge1;
/*creation of matrix from variables (n by 1 column matrix for each variable*/
read all var {'dmed_other1'} into dmed1;
read all var {'dmed_other2'} into dmed2;
read all var {'dmed_other3'} into dmed3;
read all var {'dmed_other4'} into dmed4;
read all var {'dmed_other5'} into dmed5;
read all var {'dmed_other6'} into dmed6;
read all var {'dmed_other7'} into dmed7;
read all var {'dmed_other8'} into dmed8;
read all var {'dmed_other9'} into dmed9;
read all var {'anestetics'} into anestetics;
read all var {'antibiotics'} into antibiotics;
read all var {'cardiovascular'} into cardiovascular;
read all var {'hormonal'} into hormonal;
read all var {'others'} into others;
read all var {'fluids_electrolytes'} into fluids_electrolytes;
/*matrix initiation*/
a=39; /*number of donor records*/
b=138; /*number of entries in longest list of (Anestetics, Antibiotics, Cardiovascular, Hormonal, Others)*/
c=9; /*number of open text field variable names (i.e. dmed_other 1-9) that need to be re-coded*/
/*create macro and DO loop for matrix calculations*/
%macro evaluation;
dmedrecode=j(a,c,0);
%do m=1 %to 9; /*number of open text field names that need to be re-coded*/
do n=1 to a; /*number of donors being examined for each variable*/
do p=1 to b; /*number of entries in variable list that need to be checked*/
/*coding notes: anestetics=1, antibiotics=2, cardiovascular=3, hormonal=4, others=5, fluids_electrolytes=6*/
if (dmed&m[n]=anestetics[p]) then dmedrecode[n,&m]=1;
if (dmed&m[n]=antibiotics[p]) then dmedrecode[n,&m]=2;
if (dmed&m[n]=cardiovascular[p]) then dmedrecode[n,&m]=3;
if (dmed&m[n]=hormonal[p]) then dmedrecode[n,&m]=4;
if (dmed&m[n]=others[p]) then dmedrecode[n,&m]=5;
if (dmed&m[n]=fluids_electrolytes[p]) then dmedrecode[n,&m]=6;
end;
end;
%end;
/*create matrix*/
create textmedrecode from dmedrecode;
append from dmedrecode;
close textmedrecode;
%mend evaluation;
%evaluation;
quit;
/*check to make sure all the text meds were re-coded into a category*/
proc freq data=textmedrecode;
tables col1 col2 col3 col4 col5 col6 col7 col8 col9;
run;
/*open up flags, see which meds are flagged,add to categories, then re-run*/
data flags;
merge opentextmerge1 textmedrecode;
run;
data flags1;
set flags;
if dmed_other1^="none" and col1=0 then flags1=1;
if dmed_other2^="none" and col1=0 then flags2=1;
if dmed_other3^="none" and col1=0 then flags3=1;
if dmed_other4^="none" and col1=0 then flags4=1;
if dmed_other5^="none" and col1=0 then flags5=1;
if dmed_other6^="none" and col1=0 then flags6=1;
if dmed_other7^="none" and col1=0 then flags7=1;
if dmed_other8^="none" and col1=0 then flags8=1;
if dmed_other9^="none" and col1=0 then flags9=1;
run;
/********step 6: join coded medication list into main study data table*/
DATA study5;
merge study4 textmedrecode;
run;
Out[73]:
Now we need to re-code the pre-coded donor management medications (eg steroids, diuretics, etc) and integrate them into the categories
In [75]:
DATA study6;
set study5;
pt_steroids1=0 ;
pt_diuretics1=0 ;
pt_t31=0 ;
pt_t41=0 ;
anticonv1=0 ;
antihype1=0 ;
vasodi1l=0 ;
ddavp1=0 ;
heparin1=0 ;
arginine1=0 ;
insulin1=0 ;
if pt_steroids='Y' then pt_steroids1=4 ;
if pt_diuretics='Y' then pt_diuretics1=3 ;
if pt_t3='Y' then pt_t31=4 ;
if pt_t4='Y' then pt_t41=4 ;
if anticonv='Y' then anticonv1=5 ;
if antihype='Y' then antihype1=3 ;
if vasodil='Y' then vasodi1l=3 ;
if ddavp='Y' then ddavp1=3 ;
if heparin='Y' then heparin1=5 ;
if arginine='Y' then arginine1=3 ;
if insulin='Y' then insulin1=4 ;
run;
/*donormedsreceived=0 means the donor never received a single drug during the hospital stay,
=1 means that donor received 1 or more meds*/
DATA study7;
set study6;
donormedsreceived=1;
if (pt_diuretics1=0 and pt_t31=0 and pt_t41=0 and anticonv1=0 and antihype1=0 and vasodi1l=0 and ddavp1=0 and heparin1=0 and arginine1=0 and insulin1=0 and col1=0 and col2=0 and col3=0 and col4=0 and col5=0 and col6=0 and col7=0 and col8=0 and col9=0) then donormedsreceived=0;
anestetics=0;
antibiotics=0;
cardiovascular=0;
hormonal=0;
others=0;
fluids_electrolytes=0;
if (pt_diuretics1=1 or pt_t31=1 or pt_t41=1 or anticonv1=1 or antihype1=1 or vasodi1l=1 or ddavp1=1 or heparin1=1 or arginine1=1 or insulin1=1 or col1=1 or col2=1 or col3=1 or col4=1 or col5=1 or col6=1 or col7=1 or col8=1 or col9=1 ) then anestetics=1;
if (pt_diuretics1=2 or pt_t31=2 or pt_t41=2 or anticonv1=2 or antihype1=2 or vasodi1l=2 or ddavp1=2 or heparin1=2 or arginine1=2 or insulin1=2 or col1=2 or col2=2 or col3=2 or col4=2 or col5=2 or col6=2 or col7=2 or col8=2 or col9=2 ) then antibiotics=2;
if (pt_diuretics1=3 or pt_t31=3 or pt_t41=3 or anticonv1=3 or antihype1=3 or vasodi1l=3 or ddavp1=3 or heparin1=3 or arginine1=3 or insulin1=3 or col1=3 or col2=3 or col3=3 or col4=3 or col5=3 or col6=3 or col7=3 or col8=3 or col9=3 ) then cardiovascular=3;
if (pt_diuretics1=4 or pt_t31=4 or pt_t41=4 or anticonv1=4 or antihype1=4 or vasodi1l=4 or ddavp1=4 or heparin1=4 or arginine1=4 or insulin1=4 or col1=4 or col2=4 or col3=4 or col4=4 or col5=4 or col6=4 or col7=4 or col8=4 or col9=4 ) then hormonal=4;
if (pt_diuretics1=5 or pt_t31=5 or pt_t41=5 or anticonv1=5 or antihype1=5 or vasodi1l=5 or ddavp1=5 or heparin1=5 or arginine1=5 or insulin1=5 or col1=5 or col2=5 or col3=5 or col4=5 or col5=5 or col6=5 or col7=5 or col8=5 or col9=5 ) then others=5;
if (pt_diuretics1=6 or pt_t31=6 or pt_t41=6 or anticonv1=6 or antihype1=6 or vasodi1l=6 or ddavp1=6 or heparin1=6 or arginine1=6 or insulin1=6 or col1=6 or col2=6 or col3=6 or col4=6 or col5=6 or col6=6 or col7=6 or col8=6 or col9=6 ) then fluids_electrolytes=6;
run;
PROC freq data=study7;
tables donormedsreceived anestetics antibiotics cardiovascular hormonal others fluids_electrolytes;
run;
PROC freq data=study7;
tables hospbin*donormedsreceived hospbin*anestetics hospbin*antibiotics hospbin*cardiovascular hospbin*hormonal hospbin*others hospbin*fluids_electrolytes;
run;
Out[75]:
Descriptive statistics¶
NOTE: See section 'Verification of Results' above for statistical testing and p-values for age, sex, BMI, and ethnicity.
In [44]:
PROC univariate data=study4;
var oppc_ageR
oppc_BMIR
cpeptide
transport_duration_hours;
run;
proc sort data=study4;
by hospbin;
run;
proc means data=study4 n mean std;
var oppc_ageR
oppc_BMIR
cpeptide
transport_duration_hours;
by hospbin;
run;
proc freq data=study4;
tables oppc_gender*hospbin
oppc_ethnicity*hospbin
cod_final*hospbin;
run;
Out[44]:
Need to obtain the range of hospitalization stay
In [45]:
PROC means data=study4 n mean std median min max;
var Hospitalization_Stay_Days;
run;
Out[45]:
In [46]:
proc sort data=study4;
by oppc_ageR;
run;
proc print data=study4;
run;
proc freq data=study4;
tables age_group;
run;
Out[46]: