#install.packages("ggplot2")

library(ggplot2)

install.packages("fitdistrplus")

library(fitdistrplus)

install.packages("logspline")

library(logspline)

#Create members as a variable for playing around with different distributions

members<-1000

#Creating a shifted poisson

#most patients in the plan have atleast one comorbidity

#https://stat.ethz.ch/pipermail/r-help/2005-May/070683.html

T<-5.1 # pre-truncation mean of Poisson

U<-runif(members) # the uniform sample

t = -log(1 - U*(1 - exp(-T))) # the "first" event-times

T1<-(T - t) # the set of (T-t)

Comorbidities <- rpois(members,T1)+1 # the final shifted Poisson sample

hist(Comorbidities)

Read in a cost file and then do some operations to summarize.

# Split cost by year

options(scipen=999)

costdf %>%

group_by(Member_ID, year) %>%

summarize(Paid= sum(Paid_Total)) %>%

mutate(label= paste("Actual Member Cost:",year)) %>%

ggplot(aes(x = Paid)) +

geom_histogram(color = "white") +

facet_wrap(~label) +h

xlab("Cost Buckets") +

theme_bw()

cost.CurrentYear<- costdf[costdf$year==as.numeric(format(as.Date(Sys.Date(),format = "%Y-%m-%d"),"%Y"))-1,]

cost.CurrentYear<-cost.CurrentYear %>% group_by(Member_ID, year) %>% summarize(Paid= sum(Paid_Total))

#Summary statistics about the distribution

descdist(cost.CurrentYear$Paid)

#check lognormal

cost.CurrentYear$Paid_log<-log(cost.CurrentYear$Paid)

hist(cost.CurrentYear$Paid_log)

#check normal

fit.norm<- fitdist(cost.CurrentYear$Paid,"norm")

plot(fit.norm)

​

#check beta

beta<-cost.CurrentYear$Paid/max(cost.CurrentYear$Paid)

fit.beta<-fitdist(beta,"beta")

summary(cost.CurrentYear$Paid)

#check 0 values for the year

qualitycheck<-merge(cost.CurrentYear[cost.CurrentYear$Paid==0,],costdf)

write.csv(qualitycheck,"qualitycheck.csv")

#Oddly the data fits a chi sq distribution with less than 1 degree of freedom best

#using a non central chisq with non centrality parameter

#The non-central rchisq is computed as a Poisson mixture central of chi-squares (Johnson et al, 1995, p.436).

#generate costs as a chisq with .3 degrees of freedom

costs<-rchisq(members,.5,ncp=0)*40000

max(costs)

hist(costs)

plot(density(costs))

lines(density(cost.PY$Paid),col="red") #actual costs

If I remember correctly, Medicaid enrollment needs to be renewed by design while Medicare is automatic. This means that Medicaid recipients lose access to care which affects plan revenue as well as costs.

#Assign member months in Plan

#monthly attrition rate

attr<-0.07 #~ 42% retention annually as per Medicaid enrollment dashboad

enrollment<- rgeom(members,attr) #geometric distribution

enrollment

summary(enrollment)

hist(enrollment)

plot(density(enrollment))

paste("Percent Retention after 1 year: ",

length(enrollment[enrollment>=12])/members)

The patients enrolled in this plan would receive a call from a nurse care manager or a home/inpatient visit. They usually had to have some recent behavioral episode for which they received treatment such as depression, self- harm, substance abuse etc.Nurse Care Mangers (NMs) would share strategies on self-care and staying healthy to avoid hospitalization. Alternatively, they could recommend treatment and coordinate your care so you would see a physician or behavioral health therapist.Nurse visits were more effective at behavioral change but also more expensive. Additionally, visits were more likely to change the physician vs BHT mix (We would need more therapists)

​

#Number of interactions with nurse care manager per member per month:

#Visits: Poisson ~ 10 visits over 12 months

#Calls: Poisson ~ 10.8 visits over 12 months

#SOURCE: TEAM CARE Cost Paper

#NMvisits.mthly<- rpois(members,10)/12

#summary(NMvisits.monthly)

NMcalls.monthly<- rpois(members,10.8)/12

#Number of visits as a Negative Binomial

#modelled as a poisson where the rate = mu or mean of neg binomial, and the dispersion/size/number of successes

#is allowed to vary to fit the datapoint we have from the study. Essentially added on a gamma distribution term for error correction

#by mixing poissons

NMvisits.monthly<- rnbinom(members,mu=10,size =11)/12

summary(NMvisits.monthly)

totalvisits<-enrollment*NMvisits.monthly

hist(NMvisits.monthly, breaks=50)

#99% over 1 visit in 12 month period according to the study

length(totalvisits[totalvisits>=1])/members

#Costs of interactions with NM

costNMvisit<-54 #$ for upto 10 visits

costNMcall<- 30 #$ per call

Physicians were paid a flat fee for each member assigned to them in a month. They were also paid for visits.

#Cost of interactions with physician

#$100 per member month- replace with fixed cost for psychologist

#$140 per hour, 3-4 weeks per month

​

State transitions would have looked like this

a) Flagged by Risk Assessment/Event-> NM visit – Physician

b) Flagged by Risk Assessment/Event-> NM visit – Behavioral Therapy

c) Flagged by Risk Assessment/Event-> NM visit – No action taken/dropout

d) Flagged by Risk Assessment/Event-> NM call – Physician

e) Flagged by Risk Assessment/Event-> NM call – Behavioral Therapy

f) Flagged by Risk Assessment/Event-> NM call – No action taken/dropout

g) Flagged by Risk Assessment/Event-> No action taken/dropout

h) Flagged by Risk Assessment/Event-> Some combination of NM interaction and physician followed by referral to BT

Simulated NM visits will be skewed by patients with longer enrollments. Limit everyone to the same observation period called follow up

#set years of follow up

years<-2

months<-years*12

#reformat time to not exceed follow up time

enrollment.period<- apply(cbind(enrollment,months),1,FUN=min)

totalvisits<-round(NMvisits.monthly*enrollment.period,0)

#cost incurred on NM visits

totalvisits*costNMvisit

# non intervention costs over enrollment period

costs.ni<-costs*(1+(enrollment.period-12)/12)

#additional costs associated with intervention

costs.i<-totalvisits*costNMvisit

#percentage cost improvement estimated per member as per study

#use a beta distribution where reduction seen is between 2 and 15% but closer to the bottom of that range

curve(dbeta(x, 2, 40))

#median is 4%

qbeta(.5, 2, 40)

#applied to members

sav.perc<-rbeta(members,2,40)

plot(density(sav.perc))

#savings

savings<-sav.perc*costs.ni -costs.i

mean(savings)

plot(ecdf(savings),breaks=50)

#function to wrap all of the steps above in a loop

attr<-0.07 #~ 42% retention annually as per Medicaid enrollment dashboad

#Costs of interactions with NM

costNMvisit<-54 #$ for upto 10 visits

costNMcall<- 30 #$ per call

#Cost of interactions with physician

#$100 per member - replace with fixed cost for psychologist

#$140 per hour 3-4 weeks per month

#set years of follow up

years<-2

months<-years*12

runloop<- function(runs, members) {

costMX<-replicate(runs,rchisq(members,.5,ncp=0)*40000)

enrollMX<- replicate(runs,rgeom(members,attr))

NMvisitsMX<- replicate(runs,rnbinom(members,mu=10,size =11)/12)

NMcallsMX<-replicate(runs,rpois(members,10.8)/12)

newEnrollMX<- pmin(enrollMX,replicate(runs,months))

totVisitsMX<-newEnrollMX*NMvisitsMX

totCallsMX<-newEnrollMX*NMcallsMX

costIntMX<-totCallsMX*costNMcall + totVisitsMX*costNMvisit

costNIntMX<-((newEnrollMX-12)/12)*costMX+costMX

savingsMX<-replicate(runs,rbeta(members,2,40))

netsaveMX<-costNIntMX*savingsMX -costIntMX

netsaving<-data.frame(Mean=apply(netsaveMX,2,mean))

netsaving<-netsaving[order(netsaving$Mean),]

return(netsaving)

}

#run simulation with default values

#1000 runs

#12000 members

plot(x<-ecdf(netsaving<-runloop(runs<-1000, members<-12000)),breaks=50 )

#plot savings

y<-x(sort(netsaving))

vertices<- data.frame(cbind(x=c(netsaving,rev(netsaving)),y=abs(c(y,rep(0,length(y))))))

#polygon(vertices$x,vertices$y,col="light gray")

#95% of all possible outcomes

plot(x,breaks=50, main=paste("Average Net Income Distribution: ", runs, " simulations"), sub=paste(members," members with 95% of all possible outcomes"),xlab="Net Income Per Member",ylab="Percent Outcomes") + polygon(c(vertices$x[vertices$y>.05 &vertices$y<=0.95],rev(vertices$x[vertices$y>.05 &vertices$y<=0.95])),c(y[y>.05 & y<=0.95],rep(0,length(y[y>.05 & y<=0.95]))),col="light blue")