library("rpartOrdinal") data("lowbwt") lowbwt$Category <- factor(ifelse(lowbwt$bwt<=2500,3,ifelse(lowbwt$bwt<=3000,2, ifelse(lowbwt$bwt<=3500,1,0))),ordered=TRUE) otwoing.rpart <- rpart(Category~age+lwt+race+smoke+ptl+ht+ui+ftv, data=lowbwt,method=twoing) plot(otwoing.rpart) text(otwoing.rpart,pretty=TRUE) post(otwoing.rpart,filename="TwoingLowbwt.ps",use.n=FALSE,title="", horizontal=FALSE) ordinal.rpart <- rpart(Category~age+lwt+race+smoke+ptl+ht+ui+ftv, data=lowbwt,method=ordinal) plot(ordinal.rpart) text(ordinal.rpart,pretty=TRUE) # Fit using the linear loss function linear.loss.rpart <- rpart(Category~age+lwt+race+smoke+ptl+ht+ui+ftv, data=lowbwt,method="class", parms=list(loss=loss.matrix(method="linear",lowbwt$Category))) plot(linear.loss.rpart) text(linear.loss.rpart,pretty=TRUE) # Fit using the quadratic loss function quad.loss.rpart <- rpart(Category~age+lwt+race+smoke+ptl+ht+ui+ftv, data=lowbwt,method="class", parms=list(loss=loss.matrix(method="quadratic",lowbwt$Category))) plot(quad.loss.rpart) text(quad.loss.rpart,pretty=TRUE) library(rpartOrdinal) job.satis <- factor(c(1,rep(2,3),rep(3,10),rep(4,6),rep(1,2),rep(2,3), rep(3,10),rep(4,7),1,rep(2,6),rep(3,14),rep(4,12),2,rep(3,9),rep(4,11)), ordered=TRUE,labels=c("Very Dissatisfied","Little Dissatisfied", "Moderately Satisfied","Satisfied")) income <- factor(c(rep(1,20),rep(2,22),rep(3,33),rep(4,21)),ordered=TRUE, labels=c("<15,000","15,000-25,000","25,000-40,000",">40,000")) table(job.satis,income) ordinal.gamma(job.satis,income) # V-fold cross-validation V <- 5 n <- length(lowbwt$Category) leave.out <- trunc(n/V) # sample the observations o <- sample(1:n) # the groups list stores the sampled observation numbers # belonging to each of the V folds groups <- vector("list", V) for (j in 1:(V - 1)) { jj <- (1 + (j - 1) * leave.out) groups[[j]] <- (o[jj:(jj + leave.out - 1)]) } groups[[V]] <- o[(1 + (V - 1) * leave.out):n] # Initialize fitted class vector for each method linear.fit <- rep(NA, n) quad.fit <- rep(NA, n) ordinal.fit <- rep(NA, n) twoing.fit <- rep(NA, n) for (j in 1:V) { # Fit ordinal CT to training data for Vth fold ordinal.rpart <- rpart(Category~age+lwt+race+smoke+ptl+ht+ui+ftv, data=lowbwt,subset= -groups[[j]],method=ordinal) # Get ordinal predicted class for left out fold ordinal.fit[groups[[j]]] <- predict(ordinal.rpart,newdata=lowbwt[groups[[j]],]) # Fit ordinal twoing CT to training data for Vth fold twoing.rpart <- rpart(Category~age+lwt+race+smoke+ptl+ht+ui+ftv, data=lowbwt,subset= -groups[[j]],method=twoing) # Get ordinal twoing predicted class for left out fold twoing.fit[groups[[j]]] <- predict(twoing.rpart,newdata=lowbwt[groups[[j]],]) # Fit CT with linear loss to training data for Vth fold linear.rpart <- rpart(Category~age+lwt+race+smoke+ptl+ht+ui+ftv, data=lowbwt,subset= -groups[[j]], parms=list(loss=loss.matrix(method="linear",lowbwt$Category))) # Get predicted class when using linear loss for left out fold phat <- predict(linear.rpart, newdata=lowbwt[groups[[j]],]) linear.fit[groups[[j]]] <- apply(phat,1,which.max) # Fit CT with quadratic loss to training data for Vth fold quadratic.rpart <- rpart(Category~age+lwt+race+smoke+ptl+ht+ui+ftv, data=lowbwt,subset= -groups[[j]], parms=list(loss=loss.matrix(method="quad",lowbwt$Category))) # Get predicted class when using quadratic loss for left out fold phat <- predict(quadratic.rpart, newdata=lowbwt[groups[[j]],]) quad.fit[groups[[j]]] <- apply(phat,1,which.max) } # V-fold cross-validation estimate of gamma for each method ordinal.gamma(lowbwt$Category,twoing.fit) ordinal.gamma(lowbwt$Category,ordinal.fit) ordinal.gamma(lowbwt$Category,linear.fit) ordinal.gamma(lowbwt$Category,quad.fit) library("rpartOrdinal") library("ALL") data("ALL") # Keep only patients for which BT is B1, B2, B3, or B4. BALL <- ALL[,is.element(pData(ALL)$BT,c("B1","B2","B3","B4"))] stage <- factor(pData(BALL)$BT,levels=c("B1","B2","B3","B4"),ordered=TRUE) Bcell <- data.frame(t(exprs(BALL)),stage) # Ordered twoing CT otwoing.rpart <- rpart(stage~.,data=Bcell,method=twoing) plot(otwoing.rpart) text(otwoing.rpart) # Ordinal CT ord.rpart <- rpart(stage~.,data=Bcell,method=ordinal) plot(ord.rpart) text(ord.rpart) # Generalized Gini with linear loss linear.loss.rpart <- rpart(stage~.,data=Bcell, parms=list(loss=loss.matrix(method="linear",Bcell$stage))) plot(linear.loss.rpart); ext(linear.loss.rpart) # Generalized Gini with quadratic loss quad.loss.rpart <- rpart(stage~.,data=Bcell, parms=list(loss=loss.matrix(method="quad",Bcell$stage))) plot(quad.loss.rpart) text(quad.loss.rpart) # V-fold cross-validation V <- 5 n <- length(Bcell$stage) leave.out <- trunc(n/V) # Set the random seed so that results may be replicated set.seed(123) # sample the observations o <- sample(1:n) # the groups list stores the sampled observation numbers # belonging to each of the V folds groups <- vector("list", V) for (j in 1:(V - 1)) { jj <- (1 + (j - 1) * leave.out) groups[[j]] <- (o[jj:(jj + leave.out - 1)]) } groups[[V]] <- o[(1 + (V - 1) * leave.out):n] # Initialize fitted class vector for each method linear.fit <- rep(NA, n) quad.fit <- rep(NA, n) ordinal.fit <- rep(NA, n) twoing.fit <- rep(NA, n) for (j in 1:V) { # Create training data for Vth fold train <- Bcell[-groups[[j]],] # Fit the four methods to training data ordinal.rpart <- rpart(stage~.,data=train,method=ordinal) twoing.rpart <- rpart(stage~.,data=train,method=twoing) linear.rpart <- rpart(stage~.,data=train, parms=list(loss=loss.matrix(method="linear",train$stage))) quad.rpart <- rpart(stage~.,data=train, parms=list(loss=loss.matrix(method="quad",train$stage))) rm(train) # Create test data for Vth fold test <- Bcell[groups[[j]],] # Predict class for each of the four CT methods ordinal.fit[groups[[j]]] <- predict(ordinal.rpart,newdata=test) twoing.fit[groups[[j]]] <- predict(twoing.rpart,newdata=test) phat <- predict(linear.rpart, newdata=test) linear.fit[groups[[j]]] <- apply(phat,1,which.max) rm(phat) phat <- predict(quad.rpart, newdata=test) quad.fit[groups[[j]]] <- apply(phat,1,which.max) rm(ordinal.rpart,twoing.rpart,linear.rpart,quad.rpart,phat,test) } # Five-fold CV gamma estimate ordinal.gamma(Bcell$stage,ordinal.fit) ordinal.gamma(Bcell$stage,twoing.fit) ordinal.gamma(Bcell$stage,linear.fit) ordinal.gamma(Bcell$stage,quad.fit)