Demo for identity fragmentation bias

multi_device_sims.R

# Shiny App: Demo of Identity Fragmentation Bias
# Tesary Lin and Sanjog Misra
# Dec. 7, 2021

library(shiny)
library(matrixStats)
library(repr)
library(plotrix)

# Fixed Parameters not in UI:
N = 1000 # Number of Consumers
K = 1 # Number of Covariates 
J = 2 # Number of Fragments 
# subs = FALSE # If TRUE, ignore user-specified bprob & simulate estimates under device substitution
theta = c(5,1) # model parameters dim=[1(intercept)+K(slopes)]
NR=1000 # Number of Monte Carlo reps

# Define UI for application that draws a histogram
ui <- fluidPage(
    
    # Application title
    titlePanel(""), #Fragmentation Bias: Demo
    
    sidebarLayout(
        sidebarPanel(
            numericInput("xprob","Prob. of ad exposure on desktop:", min = 0,max = 1,value = 0.5),
            numericInput("bprob","Prob. of buying on desktop:", min = 0,max = 1,value = 0.5),
            numericInput("cor","Correlation between ad exposures across devices:", min = 0, max = 1,value = 0, step = 0.05),
            radioButtons("subs", "Cross-device substitution", choices = c("FALSE" = "FALSE", "TRUE" = "TRUE"), 
                         selected = "FALSE", inline=TRUE),
            width = 2
        ),
        # Show a plot of the generated distribution
        mainPanel(
            plotOutput("paraPlot")           
        )
    )
)


# Define server logic required to draw a histogram
server <- function(input, output) {
    
    output$paraPlot <- renderPlot({
        # Function to run a single Simulation
        runsim = function(seed){
            set.seed(seed)
            # Generate true parameters
            theta.est = list()
            theta.est$tru = theta
            # Generate fragmented & true x (exposure)
            x=list()
            fullx = matrix(0,N,K)
            x[[1]] = matrix(rbinom(N*K,1,prob = rbinom(N*K,1, prob = input$xprob)),N,K,byrow=TRUE)
            fullx = fullx + x[[1]]
            if (input$cor == 0){
                corx = 0
            } else{
                corx = matrix(rbinom(N*K,1,prob = rbinom(N*K,1, prob = input$cor)),N,K,byrow=TRUE)
            }
            x[[2]] = corx * x[[1]] + (1 - corx) * matrix(rbinom(N*K,1,prob = rbinom(N*K,1,prob = (1 - input$xprob))),N,K,byrow=TRUE)
            fullx = fullx + x[[2]]
            # Generate true y (purchases)
            fully = cbind(1,fullx)%*%theta + rnorm(N)
            
            # Generate device purchase inclination
            if (input$subs == TRUE){
                bprob_s = c((x[[1]]+0.01)/(fullx+0.02))
                slam = 2 - matrix(rbinom(N*K,1,prob = rbinom(N*K,1,prob = bprob_s)),N,K,byrow=TRUE)
            }else{
                slam = sample(1:J,size=N,replace=TRUE, prob=c(input$bprob, (1 - input$bprob)))
            }
            
            # Estimates: true model
            theta.est$full = coef(lm(fully~fullx))
            
            # Estimates: device-specific effect model
            y = matrix(0,N,J)
            for(j in 1:J){
                y[which(slam==j),j]=fully[which(slam==j)]
                theta.est[[paste("frag-",j,sep="")]] = coef(lm(y[,j]~x[[j]]))
            }
            
            # Stack Fragments
            ystak = NULL
            xstak = NULL
            for(j in 1:J){
                ystak  = c(ystak,y[,j])
                xstak = rbind(xstak,x[[j]])
            }
            
            # Estimates: common effect model
            theta.est$frag = coef(lm(ystak~xstak))
            
            # return estimates
            do.call(rbind,theta.est)
        }
        
        # Monte Carlo Loop
        ints = NULL
        slopes = list()
        for(k in 1:K) slopes[[k]]=vector()
        for(r in 1:NR){
            res = runsim(r+2323)
            ints = cbind(ints,res[,1])
            for(k in 1:K){
                slopes[[k]] = cbind(slopes[[k]],res[,k+1])
            }
            cat(r,"\r")
        }
        
        # Plot parameter estimates
        color5  <- c("#000000", "#0072B2", "#009E73", "#E69F00","#F0E442") #colorblind friendly palette''
        par(mar=c(5.1, 4.1, 2.75, 2.1))
        
        # Compare slope (ad effect) estimates
        for(k in 1:K){
            plotvar = slopes[[k]]
            xr = c(min(plotvar)-.1,1.1*max(plotvar))
            options(repr.plot.width=12, repr.plot.height=8)
            plot(density(plotvar[2,], bw = 'SJ'), col=color5[2], axes=F, xlim=xr, xlab="Estimate",
                 main="Ad effect estimates comparison", lwd=3, cex.axis=1.1, cex.main = 1.75, cex.lab=1.2) #slope effect
            # python despine style axes
            box(bty="l", col="gray")
            axis(2, col="gray")
            axis(1, col="gray")
            lines(density(plotvar[3,]),col=color5[4],lwd=3)
            lines(density(plotvar[4,]),col=color5[5],lwd=3)
            lines(density(plotvar[5,]),col=color5[3],lwd=3)
            legend("top", bg="transparent", bty = "n",
                   legend=c("True value","User-level estimate","Fragmented (main)",
                            "Fragmented (by device): desktop","Fragmented (by device): mobile"),
                   col= color5, lty=c(2,rep(1,4)),lwd=2, cex=1.1)
            abline(v=theta[1+k],lty=2, lwd=2)
        }
    },
    width = 800, height = 500)
    
    
}

# Run the application 
shinyApp(ui = ui, server = server)