This paper examines network prominence in a co-prescription network as an indicator of opioid doctor shopping (i.e., fraudulent solicitation of opioids from multiple prescribers). Using longitudinal data from a large commercially insured population, we construct a network where a tie between patients is weighted by the number of shared opioid prescribers. Given prior research suggesting that doctor shopping may be a social process, we hypothesize that active doctor shoppers will occupy central structural positions in this network. We show that network prominence, operationalized using PageRank, is associated with more opioid prescriptions, higher predicted risk for dangerous morphine dosage, opioid overdose, and opioid use disorder, controlling for number of prescribers and other variables. Moreover, as a patient’s prominence increases over time, so does their risk for these outcomes, compared to their own average level of risk. Results highlight the importance of co-prescription networks in characterizing high-risk social dynamics.

To conduct social network analysis, we construct a patient co-prescription network in which a tie between patients indicates that they were prescribed one or more opioids by the same prescriber (as identified with a unique provider identification number). For this process, we omit opioid agents used exclusively or primarily for medication assisted treatment (MAT; e.g., buprenorphine). This strategy reduces concerns that network centrality measures were an artifice of network clustering due to sparsely located MAT-licensed providers in medically underserved areas. Network ties are undirected and weighted by the number of unique providers from which opioid co-prescriptions (minus MATs) were obtained. For example, if Patient A and Patient B were prescribed opioids by a set of the same three unique providers, the weight of their tie is three. Ties are also pooled across three quarters (T-2, T-1, and T) to account for potential lags in information diffusion and to offset the unnatural cut points imposed by quarterly observation.

The data also contain bipartite, or two-mode, network measures. This is accomplished using the generalized Co-HITS algorithm developed by Deng and colleagues (2014). This algorithm produces a PageRank score for patients that is based on both their own structural position and the prominence of the prescribers to which they are directly and indirectly connected. Many different weighting strategies are possible in the context of Co-HITS, but here we constrain all weights to be 1.

These are deidentified health claims from a large commercially insured population of about 19 million patients for the period of June 2015 through December 2016, covering the states of Georgia, North Carolina, Tennessee, Kentucky, Virginia, and West Virginia. Data are observed quarterly at the patient level and linked across administrative and health data. Patients are demographically representative of the US population with respect to gender and age, and representative of the commercially-insured population on all other measurable characteristics. The analysis sample contains 526,914 patients who contribute 2,107,656 quarterly observations.

*Brea L. Perry
*Stata code for PLOS ONE PONE-D-19-05832
*Aug 03, 2019

NOTE: zipcode not included in data because of HIPAA privacy regulations

use "DATA"

recode oud anyod opioidod (.=0)

*drop anyone who was never prescribed an opioid in 2016
bysort patid: egen totrx=total(numrx)
drop if totrx==0

gen logdegree=ln(degree+1)
lab var logdegree "Network degree logged"

egen stdpr=std(pagerank)
lab var stdpr "Standardized pagerank"

fre state
drop if state=="PR" | state=="-"

tab state, gen(st)

*Some patients switched states, so define state as state they lived in longest for purposes of MLM

bysort patid: egen numnc=total(st6)
bysort patid: egen numga=total(st2)
bysort patid: egen numwv=total(st13)
bysort patid: egen numky=total(st3)
bysort patid: egen numva=total(st12)
bysort patid: egen numtn=total(st11)
bysort patid: egen numoh=total(st8)

egen maxstate=rowmax(numnc-numoh)

gen appseven=.
replace appseven=1 if maxstate==numwv
replace appseven=2 if maxstate==numky
replace appseven=3 if maxstate==numva
replace appseven=4 if maxstate==numtn
replace appseven=5 if maxstate==numoh
replace appseven=6 if maxstate==numnc
replace appseven=7 if maxstate==numga
replace appseven=. if appseven==7 & numga==0

lab def appseven 1 "WV" 2 "KY" 3 "VA" 4 "TN" 5 "OH" 6 "NC" 7 "GA"
lab val appseven appseven

*Drop patients not in one of these states
keep if appseven!=.

*Generate city indicator for MLM
gen zipshort=int(zipcode/100)
egen citygrp=group(zipshort)
sum citygrp

xtset patid time

*Create percentile rank on PageRank
*bysort time: egen prpctile2=pctile(pagerank)

bysort time: egen prrank=rank(pagerank)
bysort time: egen maxrank=max(prrank)
gen prpctile=round(prrank/(maxrank/100))
drop maxrank
lab var prpctile "Percentile rank in 1% increments"

bysort time: egen biprrank=rank(birank)
bysort time: egen maxbirank=max(biprrank)
gen biprpctile=round(biprrank/(maxbirank/100))
drop maxbirank
lab var biprpctile "Percentile bipartite PageRank in 1% increments"

bysort patid: egen BPprpct=mean(prpctile)
gen WPprpct=prpctile-BPprpct
lab var BPprpct "Between Person percentile PageRank"
lab var WPprpct "Within Person percentile PageRank"

bysort patid: egen BPbiprpct=mean(biprpctile)
gen WPbiprpct=biprpctile-BPbiprpct
lab var BPbiprpct "Between Person bipartite percentile PageRank"
lab var WPbiprpct "Within Person bipartite percentile PageRank"

bysort patid: egen BPstdpr=mean(stdpr)
gen WPstdpr=stdpr-BPstdpr
lab var BPstdpr "Between Person Std PageRank"
lab var WPstdpr "Within Person Std PageRank"

bysort patid: egen BPlogdegree=mean(logdegree)
gen WPlogdegree=logdegree-BPlogdegree
lab var BPlogdegree "Between Person logged degree centrality"
lab var WPlogdegree "Within Person logged degree centrality"

gen inscat=.
replace inscat=1 if hmo==0 & pos==0
replace inscat=2 if pos==1
replace inscat=3 if hmo==1
lab def inscat 1 "Other" 2 "POS" 3 "HMO"
lab val inscat inscat
lab var inscat "Insurance, categorical"

gen maxmmerd=round(maxmme, 10)
fre maxmmerd, all
lab var maxmmerd "Max daily MME rounded to nearest ten"

gen odmme=maxmmerd
recode odmme (min/89=0)(90/max=1)
lab var odmme "MME>or equal to 90mg, high overdose potential"
fre odmme

gen divmme=maxmmerd
recode divmme (min/499=0)(500/max=1)
lab var divmme "Avg MME>or equal to 1000mg, diversion potential"
fre divmme

gen age10=age/10
lab var age10 "Age in tens of years"

gen BPprpct10= BPprpct/10
gen WPprpct10=WPprpct/10
gen BPbiprpct10=BPbiprpct/10
gen WPbiprpct10=WPbiprpct/10

gen numdocstrunc=numdocs
recode numdocstrunc (4/max=4)
lab var numdocstrunc "Truncated version of numdocs bc one supplemental model will not converge with original"

*Drop cases missing [no missing] 
egen miss=rowmiss(female age appseven cancer inscat numdocs pagerank numrx)
drop if miss>0

*Descriptives

xttab cancer
xttab odmme
xttab divmme
xttab female 
xttab lowinc 
xttab appseven
xttab inscat
xttab oud
xttab anyod
xtsum age degree pagerank bidegree birank prpctile biprpctile numdocs numrx numpharm maxmme 

fre BPprpct WPprpct
sum BPprpct WPprpct

xttab odmme if prpctile>=90
xttab odmme if prpctile>=95
xttab odmme if prpctile>=99

xttab divmme if prpctile>=90
xttab divmme if prpctile>=95
xttab divmme if prpctile>=99

*Pairwise correlations

pwcorr prpctile degree numrx numdocs numpharm avgmme, sig
oneway prpctile inscat
ttest prpctile, by(female)
ttest prpctile, by(lowinc)
ttest prpctile, by(divmme)
ttest prpctile, by(odmme)
ttest prpctile, by(anyod)
ttest prpctile, by(oud)

*Histograms

histogram numrx
histogram numpharm
histogram numdocs
histogram pagerank
histogram birank
histogram WPprpct
histogram BPprpct
histogram WPbiprpct
histogram BPbiprpct
histogram degree
histogram bidegree

*MULTIVARIATE MODELS

*FINAL MODELS
xtnbreg numrx female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.appseven, irr
margins, at(BPprpct10=(2 4 6 8 10))
margins, at(WPprpct10=(-6 -4 -2 0 2 4 6))
estat ic

xtlogit odmme female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.appseven, or
margins, at(BPprpct10=(2 4 6 8 10))
margins, at(WPprpct10=(-6 -4 -2 0 2 4 6))
estat ic

xtlogit oud female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.appseven, or
margins, at(BPprpct10=(2 4 6 8 10))
margins, at(WPprpct10=(-6 -4 -2 0 2 4 6))
estat ic

xtlogit anyod female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.appseven, or
margins, at(BPprpct10=(2 4 6 8 10))
margins, at(WPprpct10=(-6 -4 -2 0 2 4 6))
estat ic

*Sensitivity Analysis 1 - replicate after omitting patients with <= one prescriber per quarter
        *Switch to xtpoisson because no longer zero inflated

bysort patid: egen meandocs=mean(numdocs)
lab var meandocs "Mean number of unique doctors per quarter over all quarters"
 

xtpoisson numrx female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.appseven if meandocs>1, irr
estat ic

xtlogit odmme female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.appseven if meandocs>1, or
estat ic

xtlogit oud female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.appseven if meandocs>1, or
estat ic

xtlogit anyod female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.appseven if meandocs>1, or
estat ic

*Sensitivity Analysis 2 - replicate with bipartite PageRank

xtnbreg numrx female age10 cancer i.inscat numdocs BPbiprpct10 WPbiprpct10 i.appseven, irr
estat ic

*xtlogit and xtmelogit do not converge - will only converge when numdocs is right-truncated
xtlogit odmme female age10 cancer i.inscat numdocstrunc BPbiprpct10 WPbiprpct10 i.appseven, or 
estat ic

xtlogit odmme female age10 cancer i.inscat logdocs BPbiprpct10 WPbiprpct10 i.appseven, or 
estat ic

xtlogit oud female age10 cancer i.inscat numdocs BPbiprpct10 WPbiprpct10 i.appseven, or
estat ic

xtlogit anyod female age10 cancer i.inscat numdocs BPbiprpct10 WPbiprpct10 i.appseven, or
estat ic

*Sensitivity analysis 3 - with standardized raw pagerank

xtnbreg numrx female age10 cancer i.inscat numdocs BPstdpr WPstdpr i.appseven, irr
estat ic

xtlogit odmme female age10 cancer i.inscat numdocs BPstdpr WPstdpr i.appseven, or
estat ic

xtlogit oud female age10 cancer i.inscat numdocs BPstdpr WPstdpr i.appseven, or
estat ic

xtlogit anyod female age10 cancer i.inscat numdocs BPstdpr WPstdpr i.appseven, or
estat ic

*Sensitivity analysis 4 - Replicate with logged degree centrality

xtnbreg numrx female age10 cancer i.inscat numdocs BPlogdegree WPlogdegree i.appseven, irr
estat ic

xtlogit odmme female age10 cancer i.inscat numdocs BPlogdegree WPlogdegree i.appseven, or
estat ic

xtlogit oud female age10 cancer i.inscat numdocs BPlogdegree WPlogdegree i.appseven, or
estat ic

xtlogit anyod female age10 cancer i.inscat numdocs BPlogdegree WPlogdegree i.appseven, or
estat ic

*Sensitivity Analysis 5 - replicate with city-level fixed effects - FOOTNOTE ONLY (no table)

xtnbreg numrx female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.citygrp, irr
estat ic

xtlogit odmme female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.citygrp, or
estat ic

xtlogit oud female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.citygrp, or
estat ic

xtlogit anyod female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.citygrp, or
estat ic

*Reviewer comment regarding diversion potential

xtlogit divmme female age10 cancer i.inscat numdocs BPprpct10 WPprpct10 i.citygrp, or
estat ic

log close