diff --git a/Assignment7.Rmd b/Assignment7.Rmd index 105cbdf..98086b1 100644 --- a/Assignment7.Rmd +++ b/Assignment7.Rmd @@ -5,71 +5,202 @@ date: "11/30/2016" output: html_document --- -In the following assignment you will be looking at data from an one level of an online geography tutoring system used by 5th grade students. The game involves a pre-test of geography knowledge (pre.test), a series of assignments for which you have the average score (av.assignment.score), the number of messages sent by each student to other students about the assignments (messages), the number of forum posts students posted asking questions about the assignment (forum.posts), a post test at the end of the level (post.test) and whether or not the system allowed the students to go on to the next level (level.up). +In the following assignment you will be looking at data from an one level of an online geography tutoring system used by 5th grade students. +1.The game involves a pre-test of geography knowledge (pre.test) +2.A series of assignments for which you have the average score (av.assignment.score) +3.The number of messages sent by each student to other students about the assignments (messages) +4.The number of forum posts students posted asking questions about the assignment (forum.posts) +5.A post test at the end of the level (post.test) +6.Whether or not the system allowed the students to go on to the next level (level.up) ## Part I -#Upload data +#Get ready with the pacakages ```{r} +library(dplyr) +library(tidyr) +library(ggplot2) +library(rpart) +library(corrplot) +``` +#Upload data +```{r} +df<-read.csv("online.data.csv",stringsAsFactors =FALSE, header = TRUE ) ``` #Visualization ```{r} #Start by creating histograms of the distributions for all variables (#HINT: look up "facet" in the ggplot documentation) +df1<-gather(df,Vars,Value,-id,-level.up) -#Then visualize the relationships between variables +#Distribution for all the continuous variable +ggplot(df1,aes(x=Value))+ + geom_histogram()+ + facet_grid(.~Vars,scales = "free") -#Try to capture an intution about the data and the relationships +#Plot for all the discreet variable +df1$level_up<-"Level_up" +ggplot(df1,aes(x=level_up,fill=factor(level.up)))+ + geom_bar(position = "stack",width=0.2)+ + scale_color_discrete("Level_up_option") + +#Then visualize the relationships between variables +COR<-cor(df[,-c(1,7)]) + +corrplot(COR, order="AOE", method="circle", tl.pos="ld", type="lower", +tl.col="blue", tl.cex=0.8, tl.srt=20, + addCoef.col="black", addCoefasPercent = TRUE, + sig.level=0.50, insig = "blank") + +#Try to capture an intuition about the data and the relationships +#1. All the numeric variables are positively correlated +#2. post.test.score and message(student sent to others to ask questions) are magnitudinally correlated with r scoring .94 +#3. post.test.score and average assignment score are moderately to highly correlated as well with r scoring .76 +#4. message(student sent to others to ask questions) and average assignment score are moderately to highly correlated with r scoring .71 +#5. Number of posts in forum shows weak correaltion with other variables ``` #Classification tree ```{r} #Create a classification tree that predicts whether a student "levels up" in the online course using three variables of your choice (As we did last time, set all controls to their minimums) +#Try to use two model to see which one is more accurate + +#First one with these three variables:"post.test.score","messages","av.assignment.score" based on the correlation plot +#Second one with these three variables: "pre.test.score","post.test.score","messages" based on common sense + +#pre_prune<-rpart.control(maxdepth = ,minsplit = ) + +training_sample<-sample(nrow(df),nrow(df)*0.75) +df_training<-df[training_sample,] +df_testing<-df[-training_sample,] + + +set.seed(726) +c.tree1<-rpart(level.up~post.test.score+messages+av.assignment.score,data=df_training,method="class") +c.tree2<-rpart(level.up~pre.test.score+post.test.score+messages,data=df_training,method="class") + #Plot and generate a CP table for your tree +plotcp(c.tree1) +plotcp(c.tree2) +printcp(c.tree1) +printcp(c.tree2) #Generate a probability value that represents the probability that a student levels up based your classification tree +df_testing$predicted1<-predict(c.tree1,df_testing,type="prob")[,2] +df_testing$predicted2<-predict(c.tree2,df_testing,type="prob")[,2] # in order to do draw the ROC AND AUC -D1$pred <- predict(rp, type = "prob")[,2]#Last class we used type = "class" which predicted the classification for us, this time we are using type = "prob" to see the probability that our classififcation is based on. +#Last class we used type = "class" which predicted the classification for us, this time we are using type = "prob" to see the probability that our classififcation is based on. ``` + ## Part II #Now you can generate the ROC curve for your model. You will need to install the package ROCR to do this. ```{r} library(ROCR) +library(pROC)#(I am using a different package) +#Plot the curve(ROCR) +pred.detail1 <- prediction(df_testing$predicted1,df_testing$level.up) +plot(performance(pred.detail1, "tpr", "fpr"))+abline(0, 1, lty = 2) + +pred.detail2 <- prediction(df_testing$predicted2,df_testing$level.up) +plot(performance(pred.detail2, "tpr", "fpr"))+abline(0, 1, lty = 2) + +#Plot the curve(pROC) +ROC1<-roc(df_testing$level.up,df_testing$predicted1) +ROC2<-roc(df_testing$level.up,df_testing$predicted2) + +plot(ROC1,col="red") +plot(ROC2,col="blue") + +roc(df_testing$level.up,df_testing$predicted1,plot = TRUE,legacy.axes=TRUE,col="red",xlab="False Positive Percentage",ylab="True Positive Percentage",lwd=4) -#Plot the curve -pred.detail <- prediction(D1$pred, D1$level.up) -plot(performance(pred.detail, "tpr", "fpr")) -abline(0, 1, lty = 2) +roc(df_testing$level.up,df_testing$predicted2,plot = TRUE,legacy.axes=TRUE,col="blue",xlab="False Positive Percentage",ylab="True Positive Percentage",lwd=4) -#Calculate the Area Under the Curve -unlist(slot(performance(Pred2,"auc"), "y.values"))#Unlist liberates the AUC value from the "performance" object created by ROCR +#Here remember, if legacy.axes=TRUE wsa not specified, the x-axis would be specificty but not 1-specificty(False positive rate) +#Calculate the Area Under the Curve (ROCR) +unlist(slot(performance(pred.detail1,"auc"), "y.values")) #1 +unlist(slot(performance(pred.detail2,"auc"), "y.values")) #0.8276047 +#Unlist liberates the AUC value from the "performance" object created by ROCR + +#Calculate the Area Under the Curve (pROC) +auc(ROC1) # 1 +auc(ROC2) # 0.8276 #Now repeat this process, but using the variables you did not use for the previous model and compare the plots & results of your two models. Which one do you think was the better model? Why? ``` ## Part III #Thresholds ```{r} -#Look at the ROC plot for your first model. Based on this plot choose a probability threshold that balances capturing the most correct predictions against false positives. Then generate a new variable in your data set that classifies each student according to your chosen threshold. +# Look at the ROC plot for your first model. Based on this plot choose a probability threshold that balances capturing the most correct predictions against false positives. Then generate a new variable in your data set that classifies each student according to your chosen threshold. -threshold.pred1 <- +# Since tree1 is too perfect to predict,I decide to observe tree 2's ROC plot +#Find the threshold +roc.info1<-roc(df_testing$level.up,df_testing$predicted1,plot = TRUE,legacy.axes=TRUE) +roc.info2<-roc(df_testing$level.up,df_testing$predicted2,plot = TRUE,legacy.axes=TRUE) -#Now generate three diagnostics: +roc.df1<-data.frame(tpp=roc.info1$sensitivities*100, + fpp=(1-roc.info1$specificities)*100, + thresholds=roc.info1$thresholds) +head(roc.df1) -D1$accuracy.model1 <- +roc.df2<-data.frame(tpp=roc.info2$sensitivities*100, + fpp=(1-roc.info2$specificities)*100, + thresholds=roc.info2$thresholds) +head(roc.df2) # find the optimal threshoulds:0.4078812 -D1$precision.model1 <- +# since we have already given up the first +``` -D1$recall.model1 <- +```{r} +# using the new threshold to classify +df_testing$threshold.pred2 <-ifelse(df_testing$predicted2>0.4078812,"yes","no") -#Finally, calculate Kappa for your model according to: +# Now generate three diagnostics: +# Creating a confusion_matrix +confusion_matrix<-as.matrix(table(df_testing$threshold.pred2,df_testing$level.up)) -#First generate the table of comparisons -table1 <- table(D1$level.up, D1$threshold.pred1) +#create a function to calculate Accuracy, Precision and Recall(for 2 by 2 matrix only) +diagnositics_cal<-function(type="a",matrix){ + TN<-matrix[1] + FN<-matrix[3] + FP<-matrix[2] + TP<-matrix[4] + if(type=="p"){ + out<-TP/(TP+FN) + }else if(type=="r"){ + out<-TP/(TP+FP) + }else if(type=="a"){ + out<-(TP+TN)/(TP+TN+FP+FN) + }else{ + print("you need to specify what type of calculation you are using") + } + out +} + + +df_testing$accuracy.model2<-diagnositics_cal("a",confusion_matrix) +df_testing$precision.model2 <- diagnositics_cal("p",confusion_matrix) +df_testing$recall.model2 <-diagnositics_cal("r",confusion_matrix) + +#model2 accuracy is 0.808, precision is 0.9183673, recall is 0.6923077 + + +#Finally, calculate Kappa for your model according to: +#First generate the table of comparisons +table(df_testing$threshold.pred2,df_testing$level.up) #Convert to matrix -matrix1 <- as.matrix(table1) +matrix1 <- as.matrix(table(df_testing$threshold.pred2,df_testing$level.up)) + + +#Calculate kappa +kappa(matrix1, exact = TRUE)/kappa(matrix1) + +# kappa is not possible >1 +#install.packages("psych") +library(psych) +cohen.kappa(matrix1) #Calculate kappa kappa(matrix1, exact = TRUE)/kappa(matrix1) diff --git a/Assignment7.html b/Assignment7.html new file mode 100644 index 0000000..2893d1e --- /dev/null +++ b/Assignment7.html @@ -0,0 +1,717 @@ + + + + +
+ + + + + + + + + + + +In the following assignment you will be looking at data from an one level of an online geography tutoring system used by 5th grade students. 1.The game involves a pre-test of geography knowledge (pre.test) 2.A series of assignments for which you have the average score (av.assignment.score) 3.The number of messages sent by each student to other students about the assignments (messages) 4.The number of forum posts students posted asking questions about the assignment (forum.posts) 5.A post test at the end of the level (post.test) 6.Whether or not the system allowed the students to go on to the next level (level.up)
+#Get ready with the pacakages
+library(dplyr)
+##
+## Attaching package: 'dplyr'
+## The following objects are masked from 'package:stats':
+##
+## filter, lag
+## The following objects are masked from 'package:base':
+##
+## intersect, setdiff, setequal, union
+library(tidyr)
+library(ggplot2)
+library(rpart)
+library(corrplot)
+## corrplot 0.84 loaded
+#Upload data
+df<-read.csv("online.data.csv",stringsAsFactors =FALSE, header = TRUE )
+#Visualization
+#Start by creating histograms of the distributions for all variables (#HINT: look up "facet" in the ggplot documentation)
+df1<-gather(df,Vars,Value,-id,-level.up)
+
+#Distribution for all the continuous variable
+ggplot(df1,aes(x=Value))+
+ geom_histogram()+
+ facet_grid(.~Vars,scales = "free")
+## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
+
+#Plot for all the discreet variable
+df1$level_up<-"Level_up"
+ggplot(df1,aes(x=level_up,fill=factor(level.up)))+
+ geom_bar(position = "stack",width=0.2)+
+ scale_color_discrete("Level_up_option")
+
+#Then visualize the relationships between variables
+COR<-cor(df[,-c(1,7)])
+
+corrplot(COR, order="AOE", method="circle", tl.pos="ld", type="lower",
+tl.col="blue", tl.cex=0.8, tl.srt=20,
+ addCoef.col="black", addCoefasPercent = TRUE,
+ sig.level=0.50, insig = "blank")
+
+#Try to capture an intuition about the data and the relationships
+#1. All the numeric variables are positively correlated
+#2. post.test.score and message(student sent to others to ask questions) are magnitudinally correlated with r scoring .94
+#3. post.test.score and average assignment score are moderately to highly correlated as well with r scoring .76
+#4. message(student sent to others to ask questions) and average assignment score are moderately to highly correlated with r scoring .71
+#5. Number of posts in forum shows weak correaltion with other variables
+#Classification tree
+#Create a classification tree that predicts whether a student "levels up" in the online course using three variables of your choice (As we did last time, set all controls to their minimums)
+
+#Try to use two model to see which one is more accurate
+
+#First one with these three variables:"post.test.score","messages","av.assignment.score" based on the correlation plot
+#Second one with these three variables: "pre.test.score","post.test.score","messages" based on common sense
+
+#pre_prune<-rpart.control(maxdepth = ,minsplit = )
+
+training_sample<-sample(nrow(df),nrow(df)*0.75)
+df_training<-df[training_sample,]
+df_testing<-df[-training_sample,]
+
+
+set.seed(726)
+c.tree1<-rpart(level.up~post.test.score+messages+av.assignment.score,data=df_training,method="class")
+c.tree2<-rpart(level.up~pre.test.score+post.test.score+messages,data=df_training,method="class")
+
+#Plot and generate a CP table for your tree
+plotcp(c.tree1)
+
+plotcp(c.tree2)
+
+printcp(c.tree1)
+##
+## Classification tree:
+## rpart(formula = level.up ~ post.test.score + messages + av.assignment.score,
+## data = df_training, method = "class")
+##
+## Variables actually used in tree construction:
+## [1] av.assignment.score post.test.score
+##
+## Root node error: 292/750 = 0.38933
+##
+## n= 750
+##
+## CP nsplit rel error xerror xstd
+## 1 0.931507 0 1.000000 1.000000 0.045731
+## 2 0.068493 1 0.068493 0.068493 0.015110
+## 3 0.010000 2 0.000000 0.000000 0.000000
+printcp(c.tree2)
+##
+## Classification tree:
+## rpart(formula = level.up ~ pre.test.score + post.test.score +
+## messages, data = df_training, method = "class")
+##
+## Variables actually used in tree construction:
+## [1] messages post.test.score pre.test.score
+##
+## Root node error: 292/750 = 0.38933
+##
+## n= 750
+##
+## CP nsplit rel error xerror xstd
+## 1 0.527397 0 1.00000 1.00000 0.045731
+## 2 0.017123 1 0.47260 0.47260 0.036341
+## 3 0.011986 4 0.41781 0.52397 0.037794
+## 4 0.010000 6 0.39384 0.53082 0.037976
+#Generate a probability value that represents the probability that a student levels up based your classification tree
+df_testing$predicted1<-predict(c.tree1,df_testing,type="prob")[,2]
+df_testing$predicted2<-predict(c.tree2,df_testing,type="prob")[,2] # in order to do draw the ROC AND AUC
+
+#Last class we used type = "class" which predicted the classification for us, this time we are using type = "prob" to see the probability that our classififcation is based on.
+#Now you can generate the ROC curve for your model. You will need to install the package ROCR to do this.
+library(ROCR)
+## Loading required package: gplots
+##
+## Attaching package: 'gplots'
+## The following object is masked from 'package:stats':
+##
+## lowess
+library(pROC)#(I am using a different package)
+## Type 'citation("pROC")' for a citation.
+##
+## Attaching package: 'pROC'
+## The following objects are masked from 'package:stats':
+##
+## cov, smooth, var
+#Plot the curve(ROCR)
+pred.detail1 <- prediction(df_testing$predicted1,df_testing$level.up)
+plot(performance(pred.detail1, "tpr", "fpr"))+abline(0, 1, lty = 2)
+
+## integer(0)
+pred.detail2 <- prediction(df_testing$predicted2,df_testing$level.up)
+plot(performance(pred.detail2, "tpr", "fpr"))+abline(0, 1, lty = 2)
+
+## integer(0)
+#Plot the curve(pROC)
+ROC1<-roc(df_testing$level.up,df_testing$predicted1)
+## Setting levels: control = no, case = yes
+## Setting direction: controls < cases
+ROC2<-roc(df_testing$level.up,df_testing$predicted2)
+## Setting levels: control = no, case = yes
+## Setting direction: controls < cases
+plot(ROC1,col="red")
+
+plot(ROC2,col="blue")
+
+roc(df_testing$level.up,df_testing$predicted1,plot = TRUE,legacy.axes=TRUE,col="red",xlab="False Positive Percentage",ylab="True Positive Percentage",lwd=4)
+## Setting levels: control = no, case = yes
+## Setting direction: controls < cases
+
+##
+## Call:
+## roc.default(response = df_testing$level.up, predictor = df_testing$predicted1, plot = TRUE, legacy.axes = TRUE, col = "red", xlab = "False Positive Percentage", ylab = "True Positive Percentage", lwd = 4)
+##
+## Data: df_testing$predicted1 in 142 controls (df_testing$level.up no) < 108 cases (df_testing$level.up yes).
+## Area under the curve: 1
+roc(df_testing$level.up,df_testing$predicted2,plot = TRUE,legacy.axes=TRUE,col="blue",xlab="False Positive Percentage",ylab="True Positive Percentage",lwd=4)
+## Setting levels: control = no, case = yes
+## Setting direction: controls < cases
+
+##
+## Call:
+## roc.default(response = df_testing$level.up, predictor = df_testing$predicted2, plot = TRUE, legacy.axes = TRUE, col = "blue", xlab = "False Positive Percentage", ylab = "True Positive Percentage", lwd = 4)
+##
+## Data: df_testing$predicted2 in 142 controls (df_testing$level.up no) < 108 cases (df_testing$level.up yes).
+## Area under the curve: 0.9056
+#Here remember, if legacy.axes=TRUE wsa not specified, the x-axis would be specificty but not 1-specificty(False positive rate)
+
+#Calculate the Area Under the Curve (ROCR)
+unlist(slot(performance(pred.detail1,"auc"), "y.values")) #1
+## [1] 1
+unlist(slot(performance(pred.detail2,"auc"), "y.values")) #0.8276047
+## [1] 0.9056468
+#Unlist liberates the AUC value from the "performance" object created by ROCR
+
+#Calculate the Area Under the Curve (pROC)
+auc(ROC1) # 1
+## Area under the curve: 1
+auc(ROC2) # 0.8276
+## Area under the curve: 0.9056
+#Now repeat this process, but using the variables you did not use for the previous model and compare the plots & results of your two models. Which one do you think was the better model? Why?
+#Thresholds
+# Look at the ROC plot for your first model. Based on this plot choose a probability threshold that balances capturing the most correct predictions against false positives. Then generate a new variable in your data set that classifies each student according to your chosen threshold.
+
+# Since tree1 is too perfect to predict,I decide to observe tree 2's ROC plot
+#Find the threshold
+roc.info1<-roc(df_testing$level.up,df_testing$predicted1,plot = TRUE,legacy.axes=TRUE)
+## Setting levels: control = no, case = yes
+## Setting direction: controls < cases
+
+roc.info2<-roc(df_testing$level.up,df_testing$predicted2,plot = TRUE,legacy.axes=TRUE)
+## Setting levels: control = no, case = yes
+## Setting direction: controls < cases
+
+roc.df1<-data.frame(tpp=roc.info1$sensitivities*100,
+ fpp=(1-roc.info1$specificities)*100,
+ thresholds=roc.info1$thresholds)
+head(roc.df1)
+## tpp fpp thresholds
+## 1 100 100 -Inf
+## 2 100 0 0.5
+## 3 0 0 Inf
+roc.df2<-data.frame(tpp=roc.info2$sensitivities*100,
+ fpp=(1-roc.info2$specificities)*100,
+ thresholds=roc.info2$thresholds)
+head(roc.df2) # find the optimal threshoulds:0.4078812
+## tpp fpp thresholds
+## 1 100.00000 100.00000 -Inf
+## 2 99.07407 27.46479 0.1939952
+## 3 89.81481 19.71831 0.3790491
+## 4 83.33333 18.30986 0.4753695
+## 5 78.70370 13.38028 0.6112957
+## 6 69.44444 11.97183 0.7048917
+# since we have already given up the first
+# using the new threshold to classify
+df_testing$threshold.pred2 <-ifelse(df_testing$predicted2>0.4078812,"yes","no")
+
+# Now generate three diagnostics:
+# Creating a confusion_matrix
+confusion_matrix<-as.matrix(table(df_testing$threshold.pred2,df_testing$level.up))
+
+
+#create a function to calculate Accuracy, Precision and Recall(for 2 by 2 matrix only)
+diagnositics_cal<-function(type="a",matrix){
+ TN<-matrix[1]
+ FN<-matrix[3]
+ FP<-matrix[2]
+ TP<-matrix[4]
+ if(type=="p"){
+ out<-TP/(TP+FN)
+ }else if(type=="r"){
+ out<-TP/(TP+FP)
+ }else if(type=="a"){
+ out<-(TP+TN)/(TP+TN+FP+FN)
+ }else{
+ print("you need to specify what type of calculation you are using")
+ }
+ out
+}
+
+
+df_testing$accuracy.model2<-diagnositics_cal("a",confusion_matrix)
+df_testing$precision.model2 <- diagnositics_cal("p",confusion_matrix)
+df_testing$recall.model2 <-diagnositics_cal("r",confusion_matrix)
+
+#model2 accuracy is 0.808, precision is 0.9183673, recall is 0.6923077
+
+
+#Finally, calculate Kappa for your model according to:
+#First generate the table of comparisons
+table(df_testing$threshold.pred2,df_testing$level.up)
+##
+## no yes
+## no 116 18
+## yes 26 90
+#Convert to matrix
+matrix1 <- as.matrix(table(df_testing$threshold.pred2,df_testing$level.up))
+
+
+#Calculate kappa
+kappa(matrix1, exact = TRUE)/kappa(matrix1)
+## [1] 1.043684
+# kappa is not possible >1
+#install.packages("psych")
+library(psych)
+##
+## Attaching package: 'psych'
+## The following objects are masked from 'package:ggplot2':
+##
+## %+%, alpha
+cohen.kappa(matrix1)
+## Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)
+##
+## Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries
+## lower estimate upper
+## unweighted kappa 0.55 0.64 0.74
+## weighted kappa 0.55 0.64 0.74
+##
+## Number of subjects = 250
+#Calculate kappa
+kappa(matrix1, exact = TRUE)/kappa(matrix1)
+## [1] 1.043684
+#Now choose a different threshold value and repeat these diagnostics. What conclusions can you draw about your two thresholds?
+Please submit your assignment by first “knitting” your RMarkdown document into an html file and then commit, push and pull request both the RMarkdown file and the html file.
+