Session 6 Updates

Hi all,

Session 6 is done. We covered two main ways to map perceptual data - (i) using the attribute ratings (AR) method to create p-maps and joint-space maps (JSMs), and (ii) using the overall similarity (OS) approach to create multidimensional scaling (MDS) maps.We also saw some 101 stuff on positioning, definitional terms, common positioning startegies etc. The point was to get you thinking on how the mapping process could throw insights onto positioning in general, which strategy to adopt based on what criteria etc.

OK, next, what will follow is the code and snapshots of the plots that emerge from the classwork examples I did. Again, you are strongly encouraged to replicate the classwork examples at home. Copy-paste a only a few lines of code at a time after reading the comments next to each line of code. {P.S.- the statements following a '#' are for documentation purposes only and aren't executed}.So, without further ado, let us start right away:

##########################################

1. Simple Data Visualization using biplots: USArrests example.

We use USArrests data (inbuilt R dataset) to see how it can be visualized in 2 dimensions. Just copy-paste the code below onto the R console [Hit 'enter' after the last line]. Need to install package "MASS". Don't reinstall if you have already installed it previously. A package once installed lasts forever.

rm(list = ls()) # clear workspace install.packages("MASS") # install MASS package mydata = USArrests # USArrests is an inbuilt dataset pc.cr = princomp(mydata, cor=TRUE) # princomp() is core func summary(pc.cr) # summarize the pc.cr object biplot(pc.cr) # plot the pc.cr object abline(h=0); abline(v=0) # draw horiz and vertical axes

This is what the plot should look like. Click on image for larger view.

2. Code for making Joint Space maps:

I have coded a user-defined function called JSM in R. You can use it whenever you need to make joint space maps provided just by invoking the function. All it requires to work is a perceptions table and a preference rating table. First copy-paste the entire block of code below onto your R console. Those interested in reading the code, pls copy-paste line-by-line. I have put explanations in comments ('#') for what the code is doing.

## --- Build func to run simple perceptual maps --- ## JSM = function(inp1, prefs){ #JSM() func opens # inp1 = perception matrix with row and column headers # brands in rows and attributes in columns # prefs = preferences matrix par(pty="s") # set square plotting region fit = prcomp(inp1, scale.=TRUE) # extract prin compts plot(fit$rotation[,1:2], # use only top 2 prinComps type ="n", xlim=c(-1.5,1.5), ylim=c(-1.5,1.5), # plot parms main ="Joint Space map - Home-brew on R") # plot title abline(h=0); abline(v=0) # build horiz and vert axes attribnames = colnames(inp1); brdnames = rownames(inp1) # -- insert attrib vectors as arrows -- for (i1 in 1:nrow(fit$rotation)){ arrows(0,0, x1 = fit$rotation[i1,1]*fit$sdev[1], y1 = fit$rotation[i1,2]*fit$sdev[2], col="blue", lwd=1.5); text(x = fit$rotation[i1,1]*fit$sdev[1], y = fit$rotation[i1,2]*fit$sdev[2], labels = attribnames[i1],col="blue", cex=1.1)} # --- make co-ords within (-1,1) frame --- # fit1=fitfit1$x[,1]=fit$x[,1]/apply(abs(fit$x),2,sum)[1] fit1$x[,2]=fit$x[,2]/apply(abs(fit$x),2,sum)[2] points(x=fit1$x[,1], y=fit1$x[,2], pch=19, col="red") text(x=fit1$x[,1], y=fit1$x[,2], labels=brdnames, col="black", cex=1.1) # --- add preferences to map ---# k1 = 2; #scale-down factor pref = data.matrix(prefs)# make data compatible pref1 = pref %*% fit1$x[,1:2]for (i1 in 1:nrow(pref1)){ segments(0, 0, x1 = pref1[i1,1]/k1, y1 = pref1[i1,2]/k1, col="maroon2", lwd=1.25) points(x = pref1[i1,1]/k1, y = pref1[i1,2]/k1, pch=19, col="maroon2") text(x = pref1[i1,1]/k1, y = pref1[i1,2]/k1, labels = rownames(pref)[i1], adj = c(0.5, 0.5), col ="maroon2", cex = 1.1)} # voila, we're done! #} # JSM() func ends

3. OfficeStar MEXL example done on R

Goto LMS folder 'Session 6 files'. The file 'R code officestar.txt' contains the code (which I've broken up into chunks and annotated below) and the files 'officestar data1.txt' and 'officestar pref data2.txt' contain the average perceptions or attribute table and preferences table respectively.

Step 3a: Read in the attribute table into 'mydata'.

# -- Read in Average Perceptions table -- # mydata = read.table(file.choose(), header = TRUE) mydata = t(mydata) #transposing to ease analysis mydata #view the table read # extract brand and attribute names # brdnames = rownames(mydata); attribnames = colnames(mydata)

Step 3b: Read into R the preferences table into 'prefs'.

# -- Read in preferences table -- # pref = read.table(file.choose()) dim(pref) #check table dimensions pref[1:10,] #view first 10 rows

Data reading is done. You should see the data read-in as in the figure above. We can start analysis now. Finally.

Step 3c: Run Analysis

# creating empty pref dataset pref0 = pref*0; rownames(pref0) = NULL JSM(mydata, pref0) # p-map without prefs information

The above code will generate a p-map (without the preference vectors). Should look like the image below (click for larger image):

However, to make true joint-space maps (JSMs), wherein the preference vectors are overlaid atop the p-map, run the one line code below:

JSM(mydata, pref)

That is it. That one function call executes the entire JSM sequence. The result can be seen in the image below.

Again, the JSM function is generic and can be applied to *any* dataset in the input format we just saw to make joint space maps from. Am sure you'll leverage the code for animating your project datasets. Let me or Ankit know in case any assistance is needed in this regard.

4. Session 2 survey Data on Core courses:

Lookup LMS folder 'session 6 files'. Save the data and code files to your machine. Data files are 'courses data.txt' for the raw data on perceptions and courses data prefs.txt' for the preference data with student names on it. Now let the games begin.

# read in data mydata = read.table(file.choose()) # 'courses data.txt' head(mydata) # I hard coded attribute and brand names attrib.names = c("will.recommend", "persp.change", "conceptual.value.add", "practical.relevance", "interest.sustained", "difficulty.level"); brand.names = c("GSB", "INVA", "MGTO", "SAIT")

Should you try using your project data or some other dataset, you'll need to enter the brand and attribute names for that dataset in the same order in which they appear in the dataset, separately as given above.I then wrote a simple function, titled 'pmap.inp()' to denote "p-map input", to transform the raw data into a brands-attributes average peceptions table. Note that the below code is specific to the first set of columns being the preferences data.

# construct p-map input matrices using pmap.inp() func pmap.inp = function(mydata, attrib.names, brand.names){ #> pmap.inp() func opens a1 = NULL for (i1 in 1:length(attrib.names)){ start = (i1-1)*length(brand.names)+1; stop = i1*length(brand.names); a1 = rbind(a1, apply(mydata[,start:stop], 2, mean)) } # i1 loop ends rownames(a1) = attrib.names; colnames(a1) = brand.names a1 } # pmap.inp() func ends a1 = pmap.inp(mydata, attrib.names, brand.names)

And now, we're ready to run the analysis. First the p-mapo without the prefences and then the full JSM.

# now run the JSM func on data percep = t(a1[2:nrow(a1),]); percep # prefs = mydata[, 1:length(brand.names)] prefs = read.table(file.choose(), header = TRUE) # 'courses data prefs.txt' prefs1 = prefs*0; rownames(prefs1) = NULL # null preferences doc created JSM(percep, prefs1) # for p-map sans preferences

Should produce the p-map below:

And the one-line JSM run:

JSM(percep, prefs) # for p-map with preference data

Should produce the JSM below:

Follow the rest of the HW code given to run segment-wise JSMs in the same fashion.

5. Running MDS code with Car Survey Data:

The code is in 'R code HW dataset JSMs.txt' in LMS folder 'session 6 files'. The data are in 'mds car data raw.txt'. Read them in and follow the instructions here.

# --------------------- # ### --- MDS code ---- ### # --------------------- # rm(list = ls()) # clear workspace mydata = read.table(file.choose(), header = TRUE) # 'mds car data raw.txt' dim(mydata) # view dimension of the data matrix brand.names = c("Hyundai", "Honda", "Fiat", "Ford", "Chevrolet", "Toyota", "Nissan", "TataMotors", "MarutiSuzuki")

Note that I have hard-coded the brand names into 'brand.names' If you want to use this MDS code for another dataset (for your project, say) then you'll have to likewise hard-code the brand.names in.Next, I defined a function called run.mds() that takes as input the raw data and the brand names vector, runs the analysis and outputs the MDS map. Cool, or what..

### --- build user define func run.mds --- ### run.mds = function(mydata, brand.names){ # build distance matrix # k = length(brand.names) dmat = matrix(0, k, k) for (i1 in 1:(k-1)){ a1 = grepl(brand.names[i1], colnames(mydata)); for (i2 in (i1+1):k){a2 = grepl(brand.names[i2], colnames(mydata)); # note use of Regex here a3 = a1*a2; a4 = match(1, a3); dmat[i1, i2] = mean(mydata[, a4]); dmat[i2, i1] = dmat[i1, i2] } #i2 ends } # i1 ends colnames(dmat) = brand.names; rownames(dmat) = brand.names ### --- run metric MDS --- ### d = as.dist(dmat) # Classical MDS into k dimensions # fit = cmdscale(d,eig=TRUE, k=2) # cmdscale() is core MDS func fit # view results # plot solution # x = fit$points[,1]; y = fit$points[,2]; plot(x, y, xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", xlim = c(floor(min(x)), ceiling(max(x))), ylim = c(floor(min(y)), ceiling(max(y))), type="p",pch=19, col="red"); text(x, y, labels = rownames(fit$points), cex=1.1, pos=1); abline(h=0); abline(v=0)# horiz and vertical lines drawn } # run.mds func ends

Time now to finally invoke the run.mds func and get the analysis results:

# run MDS on raw data (before segmenting) run.mds(mydata, brand.names)

The resulting MDS map looks like this:

OK, that's quite a bit now for classwork replication. Let me know if any code anywhere is not running etc due to any issues.

###################################

6. Session 6 HW:

This HW is also a group submission. You will need to co-operate with the rest of our group to get it done.

• JSM based homework:
• Collect basic demographic information about your group mates - #yrs of workex, previous industry, educational qualifications, intended major etc.
• Run individual level JSM analysis on each of your team mates (and youself) using the code below (place appropriate name in student.name = c("") in that code)
• Compare the JSMs you obtain - what salient similarities and differences do you see?
• Now, using the demographic data you have collected, speculate on which demographic characteristics are best able to explain at least some of the similarities and differences you see.
• Place (i) the 4 JSms, (ii) your list of salient similarities and differences (preferably in tabular form), (iii) the demographic profile of each group member (again, in tabular form) and (iv) the subset of demographic variables that best explain the JSMs in a PPT.
• MDS based homework:
• Construct individual level MDS maps for yourself and your group members.
• Interpret them. In particular try to see what the axes might mean.
• Interpret the clusters of similar brands (brands bunched close together) in terms of what characteristics are common among them.
• Collate your (i) MDS plots, (ii) axes interpretation on each plot and (iii) similarity cluster interpretation on each plot into a PPT
• Submit one PPT for both parts of your homework. Title slide should contain group name and member names + PGIDs. Name the slide as _session6HW.pptx
• 
  
• Note: If you don't find your name in the dataset, use a friend's observation (and demographic data) instead.

Use code below to draw individual level JSM plots:

student.name = c("Himanshu") # say, student's name is Himanshu # retain only that row in the raw data which has name 'Himanshu' mydata.test = mydata[(rownames(prefs) == student.name),] # run the pmap.inp() func to build avg perceptions table a1.test = pmap.inp(mydata.test, attrib.names, brand.names); percep.test = t(a1.test[2:nrow(a1.test),]); # introduce a small perturbation lest matrix not be of full rank percep.test = percep.test + matrix(rnorm(nrow(percep.test)*ncol(percep.test))*0.01, nrow(percep.test), ncol(percep.test)); prefs.test = prefs[(rownames(prefs) == student.name),]; prefs.test # run analysis on percep.test and prefs.test JSM(percep.test, prefs.test)

Use code below to run individual level MDS plots. Just place the apropriate student.name and run.

student.name = c("Himanshu") #change student name as reqd # retain only that row in raw data with name Himanshu mydata.test = mydata[(rownames(mydata) == student.name),]; # run analysis and save result by copy pasting onto PPT run.mds(mydata.test, brand.names)

HW deadline is 9-Nov Saturday midnight. That's it for now. Contact me with queries if any.Sudhir

Session 6 HW and other Notes

Article Update:
Found this really neat article:Different ways in which you can become a data scientist. Given that Data sciences are going to become an element of competitive advantage in the medium term, I'd say that article is of interest to all MKTR students (regardless of background).

**********************************

Update: Here's an interesting article on job trends in the coming year. 6 startup trends in 2013: bootstrapping, marketing, B2B. Thought it relevant to highlight this part of the article for our MKTR course:

5) Marketing becomes as hot as tech. “By 2017, CMOs will be spending more on IT than CIOs. Driving this massive shift is the customer data that simply did not exist a decade ago.” -Ajay Agarwal, Bain Capital Ventures

6) Service marketplaces — not individual suppliers — will become the “brand.” Just as Amazon has become a leading brand for books (versus individual publishers), consumers will look to branded marketplaces for various services, such as teaching, cleaning, or construction. -Eric Chin, general partner, Crosslink

Capital

Point? What looks arcane and abstract right now (e.g., decision trees or targeting algorithms we used) is the future. **********************************

Hi all,

Your Session 6 HW (on some targeting algorithms we covered in class) requires you to go through the short caselet 'Kirin' (PDF uploaded). This HW has 4 Qs as detailed below and will use a PPT submission format (the same as for the session 5 HW).The submission deadline is 26-Dec Wednesday midnight into a dropbox.

Some background to the Qs first. Before targeting we need to do segmentation. Q2 deals with segmentation and the interpretation of segments. However, prior to segmentation, we need to know what constructs may underlie the basis variables. Q1 deals with this aspect. You did factor analysis and segmentation already in Session 5 HW. It occurs again in this HW in Q1 and Q2 but with less emphasis. The focus is more on Q3 and Q4 where we apply the randomForest and neural net algorithms respectively.

This HW also demonstrates the importance of selecting good discriminant variables. As it turns out, the discriminant variables used here are lousy and yield remarkably low predictive accuracy rates even with such sophisticated algos. The takeaway? Methods cannot alleviate deficiencies in the data beyond a point. OK, without further ado, here we go:

Questions for Session 6 HW:

• Q1.Find what constructs may underlie the basis variables. Use factor analysis, report eigenvalue scree plot & factor loadings table. Answer the following Q sub-parts:
• (i) Which variables load less than 30% onto the factor solution? (Hint: Look for Uniqueness threshold of 1-0.30 = 0.70 or above)
• (ii) ID and label the constructs you find among the variables that do load well onto the factor solution.
• Q2. Use mclust to segment the respondents. Answer the following Q parts.
• (i) What is the optimal no. of clusters?
• (ii) Report a clusplot. What is the % of variance explained by the top 2 principal components in the cluster solution?
• (iii) ID and label the segments you find.
• Q3. Split the kirin datasets into training sample (first 212 rows) and test sample (the remaining 105 rows). Train the randomForest algorithm on the training sample. Predict test sample's segment membership. Answer the following Qs:
• (i) Record predictive accuracy in both training and test samples
• (ii) Which segment appears to havwe the highest error rate?
• (iii) List the top 3 variables that best discriminate among the segments (use Mean Decrease in Accuracy metric)
• Q4. Use the split kirin datasets to try multinomial logit with neural nets on the training and test samples. Predict test sample's segment membership.
• (i) Record predictive accuracy in both training and test samples
• (ii) Which segment appears to havwe the highest error rate?
• (iii) List the top 3 variables that best discriminate among the segments (use significance as a metric here)

Some Notes on why use R:

In case you are wondering why you are being asked to bother with R, some points to note:

• Its imperative you learn how to run MKTR analyses at least once. The reason is you won't be able to effectively lead a team of people who do what you've never done. Sure, your analytics team will run the analysis but you need to have an idea of what that entails, what to expect etc.
• Folks who have any programming experience at school or at work can probably vouch for the fact that the R language is about as straightforward as, well, English. Pls implement the R code *yourself* to get a feel of the platform. Pls help your peers out who may not be as comfortable with R.
• Those who are ambivalent or undecided about using R, I say pls give it a sincere try. Its a worthwhile investment. Learning in this course is best leveraged with a basic understanding of what a versatile analytics platform does.
• Pls share problems in the code or the data that you found, workarounds that you figured out, other packages you discovered that do the same thing better/faster etc on the blog as well. R is a community based platform and it draws its strength from a distributed user and developer base.
• Those who are determined to not touch R are free to do so. No harm, no foul. Pls borrow the plots and tables from your peers, but interpret them yourself.

Dassit for now. Gotta get to work on tomorrow's session slides and associated blog-code.

Sudhir

Targeting in R - Session 6 PDA Caselet

Update: NYT Video link for social media as a focus group

*******************************

Hi all,

This is the code for classwork MEXL example "Conglomerate's PDA". This is the roadmap for what we are going to do:

• First we segment the customer base using model based clustering or mclust, the recommended method.
• Then we randomly split the dataset into training and test samples. The test sample is about one-third of the original dataset in size, following accepted practice.
• Then we try to establish via the training sample, how the discriminant variables relate to segment membership. This is where we train the Targeting algorithm to learn about how discriminant variables relate to segment memberships.
• Then comes the real test - validate algorithm performance on the test dataset. We compare prediction accuracy across traditional and proposed methods.
• Since R is happening, there are many targeting algorithms to choose from on R. I have decided to go with one that has shown good promise of late - the randomForest algorithm. Where we had seen decision trees in Session 5, think now of 'decision forests' in a sense...
• Other available algorithms that we can run (provided there is popular demand) are artificial neural nets (multi-layer perceptrons) and Support vector machines. But for now, these are not part of this course.

So without further ado, let me start right away.

1. Segment the customer Base.

To read-in data, directly save and use the 'basis' and 'discrim' notepads I have sent you by email. Then ensure you have packages 'mclust' and 'cluster' installed before running the clustering code.

# read-in basis and discrim variables basis = read.table(file.choose(), header=TRUE) dim(basis); basis[1:3,] summary(basis) discrim = read.table(file.choose(), header=TRUE) dim(discrim); discrim[1:3,] summary(discrim) # Run segmentation on the basis.training dataset library(mclust) #invoke library fit <- Mclust(basis) # run mclust fit # view result classif = fit$classification # print cluster sizes for (i1 in 1:max(classif)){print(sum(classif==i1))} # Cluster Plot against 1st 2 principal components require(cluster) fit1=cbind(classif) rownames(fit1)=rownames(basis) clusplot(basis, fit1, color=TRUE, shade=TRUE,labels=2, lines=0)

The segmentation produces 4 optimal clusters. Below is the clusplot where, interestingly, despite our using 15 basis variables, we see decent separation among the clusters in the top 2 principal components directly.

Click on the above image for larger size.

2. Split dataset into Training & Test samples

I wrote a function that randomly splits the dataset into training and test samples. Copy-paste it as a block into the R console and we can then invoke it peacefully downstream anywhere in this R session.

# Function to split dataset into training and test samples split.data <- function(data){ #func begins train.index = sample(1:nrow(data), round((2/3)*nrow(data))) a2=1:nrow(data); test.index=a2[!(a2 %in% train.index)] ### split sample into 2/3rds training and 1/3rds test train = data[train.index,] test = data[test.index,] outp = list(train, test) outp } # func ends # run split.data func on our data data = cbind(basis, discrim, classif) outp = split.data(data) # create training and test datasets n1 = ncol(basis); n2 = ncol(data) basis.training = outp[[1]][,1:n1] discrim.training = outp[[1]][,(n1+1):(n2-1)] y.training = outp[[1]][,n2] y.test = outp[[2]][,n2] discrim.test = outp[[2]][,(n1+1):(n2-1)]

There, now we have the training and test samples defined both for basis and for discriminant matrices. Time to train the machine on the training dataset next.

3. Train the Targeting Algorithm

Before proceeding further, a quick note on what random forests are and how they relate to decision trees. Recall what decision trees are (loosely speaking) - 'classifiers' of sorts that pick the variable that best splits a given Y variable into 2 or more branches. For instance, in our car example in class, we saw that our Y variable 'Mileage' was best split first by 'Price' and then by 'Type'. We call this procedure 'recursive partitioning' and it is repeated at each node (or branch ending).

For a small set of X variables, individual decision trees are fine. But what if we had X variables numbering the dozens, scores, hundreds even? To solve this issue, instead of one decision tree, imagine we had 1000s of trees, each acting on random subsets of the original variable set. Each time a forest grows, we then summarize the fit obtained across different variables. Then, we grow another forest. Then another. That, in short, is the random forest algorithm. A forest of decision trees that allows us to efficiently parallelize some recursive partitioning tasks, keep track of variable importance and predictive fit, and thereby tackle large datasets efficiently. Again, this is just an intuitive explanation, it is not rigorous in a technical sense, so don;t take it too literally.

Ensure you have library 'randomForest' installed before trying the below code.

### --- using the randomForest algo ---- ### library(randomForest) test.rf <- randomForest(factor(y.training) ~., data = discrim.training, xtest = discrim.test, ytest = factor(y.test), ntree = 3000, proximity = TRUE, importance = TRUE) print(test.rf)

That was it. If no error showed up, the algo has executed successfully. Those few lines of code invoked a rather powerful classification algorithm that finds use in myriad scientific explorations today. Fortunately for us, its is right up our street in the Targeting arena. Let us proceed further and see how the algo has performed on in-sample and out of-sample prediction. This is the output I got on executing the above code:

The output is as straightforward to read as the input was to run. Let me summarize in bullet points for ease of interpretation.

• The output tells us that the Type of run was "Classification" implying that a regression based run for continuous data is very much possible (recall we did the regression based decision tree only for the Mileage example).
• The no. of variables tried at each split is basically how many levels our Y variable had. Well, we have 4 segments for which we want to classify customers into, so we it is 4.
• The confusion matrix is a cross-tab between the predicted and the observed classifications. Higher the numbers along the diagonal, the better. Interestingly, we have an in-sample accuracy of !62% but a test sample predictive accuracy of ~72%! Usually its the other way round.
• Does a mere 72% accuracy seem too low for such a much-touted algorithm? Aah. well, let's see. First, a random coin-toss would give us ~ 25% accuracy (roughly, across 4 classes).
• Second, there were only 107 training cases for 17 discriminant variables. If you try any parametric model like, say, Logit for this, the chances the model will give up midway are quite high indeed. And even if it goes through, don't expect much by way of variable significance in the results. Why talk, let me show you an application of the multinomial Logit for this dataset shortly.
• But what about variable significance (or 'importance') in random Forests, first? Sure, we're going there next.

4. Interpreting random Forest Output

Copy paste the following code to get variable importance tables and an associated plot. Also, the actual versus predicted segment membership tables can also be obtained using the code below.

# List the importance of the variables. rn <- round(importance(test.rf), 2) rn[order(rn[, (ncol(rn)-1)], decreasing=TRUE),] varImpPlot(test.rf) # plot variable importance # -- predictions --- # observed and predicted y outputted yhat.train = cbind(y.training, test.rf$predicted)# insample yhat.test = cbind(y.test, test.rf$test[[1]])# test sample

This is a variable importance plot. Useful to use the top few variables on the left side plot.

That was it. We used a sophisticated Targeting algorithm for our decidedly tiny and simple data application. This algo scales up massively and for many more complex applications, BTW.

5. A neural-net based Logit Model

The Q may well arise, how do I know how good the proposed new algo really is unless I run it on something known and traditional? So, how about we run it on the well-established, parametric Logistic Regression model? Since our Y has 4 levels, we would ordinarily go for the Multinomial Logit model (in the R package 'mlogit').

Let me tell you what happens when I try to replicate the above in mlogit. The program shows error messages galore. Says that for 160 rows and 17 variables, we get a matrix close to Singular (expected, BTW). So I try reducing the number of X variables. The problem persists. Finally I give up.

And turn to artificial neural-nets.

Turns out that the R package 'nnet' implements a multinomial Logit on an underlying neural net. The function multinom() is the key here. The code below will invoke, execute and help interpret multinomial Logit results for our PDA dataset.

### --- doing multinom for Targetting --- ### library(nnet) # use neural nets to fit a MNL! library(car) # for Anova func call # Build a Regression model. a1 <- multinom(formula = y.training ~ ., data = discrim.training, trace = FALSE, maxit = 1000) # For inference summary(a1, Wald.ratios=TRUE) cat('==== ANOVA ====') print(Anova(a1))

This is what I get from the ANOVA table:

Notice right off how variable significance seems to be such a problem. Well, in predictive analytics, we're more concerned with how well the variables as a group are able to predict segment membership than how significant any particular variable is. Hence, the block of code that follows gives us how well neural nets are able to predict both within sample (on the training dataset) and on the test dataset.

### Model fit in training sample ### a3 = round(a1$fitted.values,3) # probability table # build predicted Y values vector yhat = matrix(0, nrow(a3), 1) for (i1 in 1:nrow(a3)){ for (i2 in 1:ncol(a3)){ if (a3[i1, i2] == max(a3[i1,])) {yhat[i1,1] = i2}}} # confusion matrix for insample table(y.training, yhat) # predictive accuracy sum(diag(table(y.training, yhat)))/nrow(discrim.training)

This is what I got:

Wow or what? neuralnets netted us a 76% accuracy in training sample prediction. In contrast, we got only some 62% accuracy in random forests. Does that mean that neural nets are better than random forests in this application? Well, hold on.

The test dataset is still there. How well does the neural net do there? Its not necessary it will trump random forests here just because it did so in the training sample.

# build predicted Y values vector for test dataset yhat1 = predict(a1, cbind(y.test, discrim.test), type="class") # confusion matrix for multinom test dataset table(y.test, yhat1) # predictive accuracy sum(diag(table(y.test,yhat1)))/nrow(discrim.test)

See what I mean? neural nets get a 64% accuracy compared to 72% in random Forests in the test dataset. I'd say, in any given app, try both. Use whichever gives you a higher test dataset prediction accuracy.

Well, Dassit for now. See you in class tomorrow.

Sudhir