Scotch_Pref.R

# File-Name:       Scotch_Pref.R                 
# Date:            2009-11-29                                
# Author:          Drew Conway                                       
# Purpose:         Display one-dimensional item response for scotch whiskey preference
# Data Used:       whiskey, package=flexmix
# Packages Used:   zelig,ggplot
# Output File:     scotch_pref.png
# Data Output:     
# Machine:         Drew Conway's MacBook 

####                INFO                        ####
# Using an incidence matrix of scotch consumption, #
# use the one-dimension item response function in  #
# Zelig to estimate consumer preference for each   #
# brand along a single dimension. Plot these       #
# prefereces with ggplot2 such that brands         #
# estimated to be more likely consumed by the same #
# scotch drinker will be closer on the y-axis      #
        
# Load packages and data                                                            
library(Zelig)
library(ggplot2)
data(whiskey,package="flexmix")

# Create data frame from scotch preference incidence matrix
scotch.pref<-data.frame(whiskey$Incidence)

# Create one-dimensional item response with MCMC via Gibbs sampler
scotch.irt1d<-zelig(cbind(Singleton,Knockando,White.Horse,Scoresby.Rare,Ushers,Macallan,
    Grant.s,Passport,Black...White,Clan.MacGregor,Ballantine,Pinch..Haig.,
    Other.brands,Cutty.Sark,Glenfiddich,Glenlivet,J.B,Dewar.s.White.Label,
    Johnnie.Walker.Black.Label,Johnnie.Walker.Red.Label,Chivas.Regal)~NULL,model="irt1d",data=scotch.pref)

# Extract mean values from samples to estimate scotch 1-d position and create data frame
scotch.means<-c()
for(i in 1:ncol(scotch.irt1d$coeff)) {scotch.means<-append(scotch.means,mean(scotch.irt1d$coeff[,i]))}
scotch.data<-data.frame(scotch.means)
rownames(scotch.data)<-names(scotch.pref)

# Plot brands with proximity based on consumption similarity, i.e. brands 
# estimated to be more likely consumed by the same scotch drinker will
# be closer on the y-axis
png("scotch_pref.png",width=700,height=650)
scotch.plot<-ggplot(scotch.data,aes(x=0,y=scotch.means,colour=scotch.means,
    size=8,label=rownames(scotch.data)))
scotch.plot<-scotch.plot+scale_x_continuous(limits=c(-1,1))+geom_text(position=position_jitter(width=0.35,height=0.1))
scotch.plot+xlab("Scotch Brands")+ylab("Mean Values of IRT1D for Scotch Consumption Incidence")+opts(title=
    "1D Preference Estimation for Scotch Whiskey\nBased on Consumption Incidence")
dev.off()