Computation using a Matern Kernel
Julian Faraway
[1] "Created: Wed Sep 3 16:38:19 2014"We repeat the computation except using a different kernel. We choose the Matern Kernel of order 3/2 as this is perhaps the next most popular choice after the squared exponential.
load("lcdb.rda")
source("funcs.R")
GPR <- TRUESet the kernel type and choose an invers ewidth which gives similar fits to the squared exponential kernel.
kerntype <- "matern"
wvec <- 2e-2Now compute all the measures. Only those measures that depend on the GP fit will be changed. The Richards measures are unchanged.
firstdate <- 53464
daterange <- 2764
detection.limit <- 20.5Compute all the measures
uids <- lcdbi$id
nuids <- length(uids)
statmat <- NULL
types <- NULL
ids <- NULL
for(i in 1:nuids){
olc <- lcdb[[i]]
type <- as.character(lcdbi$type[i])
id <- as.character(lcdbi$id[i])
aa <- lcstats(olc,GPR,kerntype,wvec,daterange=daterange,firstdate=firstdate,detection.limit=detection.limit)
if(is.na(aa[1])) next
statmat <- rbind(statmat,aa)
types <- c(types,type)
ids <- c(ids,id)
}
colnames(statmat) <- c("meds", "iqr","shov", "maxdiff", "dscore","wander", "moveloc","nobs", "totvar", "quadvar", "famp", "fslope","trend",
"outl", "skewres", "std", "shapwilk", "lsd", "gscore", "mdev", "gtvar", "rm.amplitude", "rm.mean",
"rm.stddev", "rm.beyond1std", "rm.fpr20", "rm.fpr35", "rm.fpr50",
"rm.fpr65", "rm.fpr80", "rm.lintrend", "rm.maxslope", "rm.mad",
"rm.medbuf", "rm.pairslope", "rm.peramp", "rm.pdfp", "rm.skew",
"rm.kurtosis", "rm.std", "rm.rcorbor")
cmdb <- data.frame(statmat,type=types,id=ids,row.names=ids)
cmdb$lsd[is.na(cmdb$lsd) | is.infinite(cmdb$lsd)] <- min(cmdb$lsd[!is.infinite(cmdb$lsd)],na.rm=TRUE)-1
cmdb$fslope[is.infinite(cmdb$fslope)] <- 1Split the data into a training(2/3) and a test(1/3) sample in the same way as before. The split will not be identical to that used on the Richards measures because of the problem that not all the Richards stats are computed on the NV group resulting in the discard of some objects from the calculation.
set.seed(123)
n <- nrow(cmdb)
isel <- sample(1:n,round(n/3))
trains <- cmdb[-isel,]
tests <- cmdb[isel,]There are 1240 observations in the test set and 2480 observations in the training set.
Redo the classification for this set of measures using the Matern kernel:
opts_chunk$set(comment=NA, warning=FALSE, message=FALSE, error=FALSE)
load("feat.rda")
source("funcs.R")
require(MASS)
require(nnet)
require(ggplot2)
require(rpart)
require(rpart.plot)
require(xtable)
require(kernlab)
require(randomForest)Set up the predictors that we will use throughout. Note that I have transformed some of the variables as seen in the setup.
predform <- "shov + maxdiff + dscore + log(totvar) + log(quadvar) + log(famp) + log(fslope) + log(outl) + gscore + lsd + nudlog(gtvar) + rm.amplitude + rm.beyond1std + rm.fpr20 +rm.fpr35 + rm.fpr50 + rm.fpr80 + log(rm.maxslope) + rm.mad +asylog(rm.medbuf) + rm.pairslope + log(rm.peramp) + log(rm.pdfp) + rm.skew + log(rm.kurtosis)+ rm.std + dublog(rm.rcorbor)"
preds <- c("shov","maxdiff","dscore","totvar","quadvar","famp","fslope","outl","gscore","lsd","gtvar","rm.amplitude","rm.beyond1std","rm.fpr20","rm.fpr35","rm.fpr50","rm.fpr80","rm.maxslope","rm.mad","rm.medbuf","rm.pairslope","rm.peramp","rm.pdfp","rm.skew","rm.kurtosis","rm.std","rm.rcorbor")
tpredform <- paste(preds,collapse="+")Creation of model formulae and subset testing and training sets
Formula for classification of all types:
allform <- as.formula(paste("type ~",predform))transients vs. non-variables:
trains$vtype <- ifelse(trains$type == "nv","nv","tr")
trains$vtype <- factor(trains$vtype)
tests$vtype <- ifelse(tests$type == "nv","nv","tr")
tests$vtype <- factor(tests$vtype)
vtform <- as.formula(paste("vtype ~",predform))transients only,
trains$ttype <- trains$type
trains$ttype[trains$ttype == "nv"] <- NA
trains$ttype <- factor(trains$ttype)
tests$ttype <- tests$type
tests$ttype[tests$ttype == "nv"] <- NA
tests$ttype <- factor(tests$ttype)
trform <- as.formula(paste("ttype ~",predform))cmat <- matrix(NA,nrow=4,ncol=5)
dimnames(cmat) <- list(c("All","TranNoTran","Tranonly","Heirarch"),c("LDA","RPart","SVM","NN","Forest"))All types
LDA
Linear Discriminant analysis using the default options.
We produce the cross-classification between predicted and observed class. Note that the default priors are the proportions found in the training set.
ldamod <- lda(allform ,data=trains)
pv <- predict(ldamod, tests)
cm <- xtabs( ~ pv$class + tests$type)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 29 | 2 | 0 | 3 | 0 | 3 | 0 | 8 |
| blazar | 0 | 24 | 13 | 13 | 0 | 2 | 0 | 8 |
| cv | 0 | 3 | 85 | 32 | 1 | 2 | 1 | 14 |
| downes | 1 | 4 | 14 | 33 | 0 | 8 | 9 | 1 |
| flare | 1 | 0 | 2 | 3 | 11 | 10 | 1 | 1 |
| nv | 12 | 2 | 13 | 29 | 11 | 525 | 0 | 17 |
| rrlyr | 0 | 2 | 1 | 20 | 0 | 1 | 92 | 0 |
| sn | 1 | 2 | 17 | 2 | 1 | 7 | 0 | 143 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 66 | 5 | 0 | 2 | 0 | 1 | 0 | 4 |
| blazar | 0 | 62 | 9 | 10 | 0 | 0 | 0 | 4 |
| cv | 0 | 8 | 59 | 24 | 4 | 0 | 1 | 7 |
| downes | 2 | 10 | 10 | 24 | 0 | 1 | 9 | 1 |
| flare | 2 | 0 | 1 | 2 | 46 | 2 | 1 | 1 |
| nv | 27 | 5 | 9 | 21 | 46 | 94 | 0 | 9 |
| rrlyr | 0 | 5 | 1 | 15 | 0 | 0 | 89 | 0 |
| sn | 2 | 5 | 12 | 1 | 4 | 1 | 0 | 74 |
The overall classification rate is 0.7597.
Recursive Partitioning
roz <- rpart(allform ,data=trains)
rpart.plot(roz,type=1,extra=1)
pv <- predict(roz,newdata=tests,type="class")
cm <- xtabs( ~ pv + tests$type)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 27 | 2 | 3 | 6 | 1 | 6 | 0 | 7 |
| blazar | 0 | 17 | 6 | 19 | 0 | 1 | 4 | 2 |
| cv | 1 | 6 | 80 | 35 | 0 | 0 | 2 | 3 |
| downes | 0 | 6 | 8 | 9 | 0 | 2 | 1 | 5 |
| flare | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| nv | 11 | 1 | 19 | 50 | 22 | 528 | 3 | 32 |
| rrlyr | 0 | 4 | 1 | 13 | 0 | 7 | 93 | 5 |
| sn | 5 | 3 | 28 | 3 | 1 | 14 | 0 | 138 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 61 | 5 | 2 | 4 | 4 | 1 | 0 | 4 |
| blazar | 0 | 44 | 4 | 14 | 0 | 0 | 4 | 1 |
| cv | 2 | 15 | 55 | 26 | 0 | 0 | 2 | 2 |
| downes | 0 | 15 | 6 | 7 | 0 | 0 | 1 | 3 |
| flare | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| nv | 25 | 3 | 13 | 37 | 92 | 95 | 3 | 17 |
| rrlyr | 0 | 10 | 1 | 10 | 0 | 1 | 90 | 3 |
| sn | 11 | 8 | 19 | 2 | 4 | 3 | 0 | 72 |
The overall classification rate is 0.7194.
Support Vector Machines
Use the default choice of setting from the kernlab R package for this:
svmod <- ksvm(allform, data=trains)Using automatic sigma estimation (sigest) for RBF or laplace kernel pv <- predict(svmod, tests)
cm <- xtabs( ~ pv + tests$type)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 29 | 3 | 0 | 4 | 0 | 1 | 0 | 1 |
| blazar | 0 | 23 | 4 | 7 | 0 | 0 | 0 | 0 |
| cv | 1 | 5 | 102 | 36 | 1 | 3 | 1 | 13 |
| downes | 1 | 2 | 8 | 43 | 0 | 5 | 8 | 0 |
| flare | 0 | 0 | 0 | 1 | 7 | 4 | 0 | 0 |
| nv | 12 | 0 | 8 | 33 | 14 | 540 | 1 | 20 |
| rrlyr | 0 | 2 | 0 | 7 | 0 | 1 | 93 | 0 |
| sn | 1 | 4 | 23 | 4 | 2 | 4 | 0 | 158 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 66 | 8 | 0 | 3 | 0 | 0 | 0 | 1 |
| blazar | 0 | 59 | 3 | 5 | 0 | 0 | 0 | 0 |
| cv | 2 | 13 | 70 | 27 | 4 | 1 | 1 | 7 |
| downes | 2 | 5 | 6 | 32 | 0 | 1 | 8 | 0 |
| flare | 0 | 0 | 0 | 1 | 29 | 1 | 0 | 0 |
| nv | 27 | 0 | 6 | 24 | 58 | 97 | 1 | 10 |
| rrlyr | 0 | 5 | 0 | 5 | 0 | 0 | 90 | 0 |
| sn | 2 | 10 | 16 | 3 | 8 | 1 | 0 | 82 |
The overall classification rate is 0.8024.
Neural Net
Use the multinom() function from the nnet R package. Might work better with some scaling.
svmod <- multinom(allform, data=trains, trace=FALSE, maxit=1000, decay=5e-4)
pv <- predict(svmod, tests)
cm <- xtabs( ~ pv + tests$type)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 32 | 2 | 0 | 3 | 0 | 1 | 0 | 3 |
| blazar | 1 | 24 | 4 | 11 | 0 | 1 | 0 | 6 |
| cv | 3 | 4 | 93 | 27 | 0 | 5 | 3 | 19 |
| downes | 0 | 6 | 18 | 59 | 0 | 11 | 7 | 1 |
| flare | 0 | 0 | 0 | 1 | 7 | 4 | 0 | 1 |
| nv | 7 | 0 | 9 | 25 | 12 | 531 | 0 | 14 |
| rrlyr | 0 | 0 | 0 | 7 | 2 | 0 | 93 | 0 |
| sn | 1 | 3 | 21 | 2 | 3 | 5 | 0 | 148 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 73 | 5 | 0 | 2 | 0 | 0 | 0 | 2 |
| blazar | 2 | 62 | 3 | 8 | 0 | 0 | 0 | 3 |
| cv | 7 | 10 | 64 | 20 | 0 | 1 | 3 | 10 |
| downes | 0 | 15 | 12 | 44 | 0 | 2 | 7 | 1 |
| flare | 0 | 0 | 0 | 1 | 29 | 1 | 0 | 1 |
| nv | 16 | 0 | 6 | 19 | 50 | 95 | 0 | 7 |
| rrlyr | 0 | 0 | 0 | 5 | 8 | 0 | 90 | 0 |
| sn | 2 | 8 | 14 | 1 | 12 | 1 | 0 | 77 |
The overall classification rate is 0.796.
Random Forest
Use the randomForest package with the default settings:
tallform <- as.formula(paste("type ~",tpredform))
fmod <- randomForest(tallform, data=trains)
pv <- predict(fmod, newdata=tests)
cm <- xtabs( ~ pv + tests$type)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 30 | 1 | 0 | 3 | 0 | 1 | 0 | 2 |
| blazar | 0 | 21 | 3 | 7 | 0 | 0 | 0 | 0 |
| cv | 3 | 7 | 92 | 27 | 0 | 3 | 1 | 15 |
| downes | 0 | 4 | 17 | 57 | 0 | 7 | 8 | 0 |
| flare | 0 | 0 | 0 | 2 | 8 | 1 | 0 | 0 |
| nv | 8 | 1 | 10 | 31 | 15 | 540 | 1 | 18 |
| rrlyr | 0 | 2 | 0 | 5 | 0 | 1 | 93 | 0 |
| sn | 3 | 3 | 23 | 3 | 1 | 5 | 0 | 157 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 68 | 3 | 0 | 2 | 0 | 0 | 0 | 1 |
| blazar | 0 | 54 | 2 | 5 | 0 | 0 | 0 | 0 |
| cv | 7 | 18 | 63 | 20 | 0 | 1 | 1 | 8 |
| downes | 0 | 10 | 12 | 42 | 0 | 1 | 8 | 0 |
| flare | 0 | 0 | 0 | 1 | 33 | 0 | 0 | 0 |
| nv | 18 | 3 | 7 | 23 | 62 | 97 | 1 | 9 |
| rrlyr | 0 | 5 | 0 | 4 | 0 | 0 | 90 | 0 |
| sn | 7 | 8 | 16 | 2 | 4 | 1 | 0 | 82 |
The overall classification rate is 0.8048.
Transients vs. non-variables
LDA
Linear Discriminant analysis using the default options.
We produce the cross-classification between predicted and observed class. Note that the default priors are the proportions found in the training set.
ldamod <- lda(vtform ,data=trains)
pv <- predict(ldamod, tests)
cm <- xtabs( ~ pv$class + tests$vtype)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| nv | tr | |
|---|---|---|
| nv | 523 | 84 |
| tr | 35 | 598 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| nv | tr | |
|---|---|---|
| nv | 94 | 12 |
| tr | 6 | 88 |
The overall classification rate is 0.904.
Recursive Partitioning
roz <- rpart(vtform ,data=trains)
rpart.plot(roz,type=1,extra=1)
pv <- predict(roz,newdata=tests,type="class")
cm <- xtabs( ~ pv + tests$vtype)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| nv | tr | |
|---|---|---|
| nv | 513 | 99 |
| tr | 45 | 583 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| nv | tr | |
|---|---|---|
| nv | 92 | 15 |
| tr | 8 | 85 |
The overall classification rate is 0.8839.
Support Vector Machines
Use the default choice of setting from the kernlab R package for this:
svmod <- ksvm(vtform, data=trains)Using automatic sigma estimation (sigest) for RBF or laplace kernel pv <- predict(svmod, tests)
cm <- xtabs( ~ pv + tests$vtype)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| nv | tr | |
|---|---|---|
| nv | 524 | 65 |
| tr | 34 | 617 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| nv | tr | |
|---|---|---|
| nv | 94 | 10 |
| tr | 6 | 90 |
The overall classification rate is 0.9202.
Neural Net
Use the multinom() function from the nnet R package. Might work better with some scaling.
svmod <- multinom(vtform, data=trains, trace=FALSE, maxit=1000, decay=5e-4)
pv <- predict(svmod, tests)
cm <- xtabs( ~ pv + tests$vtype)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| nv | tr | |
|---|---|---|
| nv | 520 | 66 |
| tr | 38 | 616 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| nv | tr | |
|---|---|---|
| nv | 93 | 10 |
| tr | 7 | 90 |
The overall classification rate is 0.9161.
Random Forest
Use the randomForest package with the default settings:
tallform <- as.formula(paste("vtype ~",tpredform))
fmod <- randomForest(tallform, data=trains)
pv <- predict(fmod, newdata=tests)
cm <- xtabs( ~ pv + tests$vtype)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| nv | tr | |
|---|---|---|
| nv | 517 | 61 |
| tr | 41 | 621 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| nv | tr | |
|---|---|---|
| nv | 93 | 9 |
| tr | 7 | 91 |
The overall classification rate is 0.9177.
Transients only
LDA
Linear Discriminant analysis using the default options.
We produce the cross-classification between predicted and observed class. Note that the default priors are the proportions found in the training set.
ldamod <- lda(trform ,data=trains)
pv <- predict(ldamod, tests)
cm <- xtabs( ~ pv$class + tests$ttype)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | rrlyr | sn | |
|---|---|---|---|---|---|---|---|
| agn | 37 | 3 | 4 | 7 | 4 | 0 | 11 |
| blazar | 0 | 23 | 7 | 13 | 0 | 0 | 6 |
| cv | 1 | 3 | 91 | 27 | 1 | 1 | 13 |
| downes | 3 | 7 | 17 | 65 | 5 | 9 | 2 |
| flare | 1 | 0 | 1 | 3 | 11 | 1 | 1 |
| rrlyr | 0 | 0 | 0 | 17 | 0 | 92 | 0 |
| sn | 2 | 3 | 25 | 3 | 3 | 0 | 159 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | rrlyr | sn | |
|---|---|---|---|---|---|---|---|
| agn | 84 | 8 | 3 | 5 | 17 | 0 | 6 |
| blazar | 0 | 59 | 5 | 10 | 0 | 0 | 3 |
| cv | 2 | 8 | 63 | 20 | 4 | 1 | 7 |
| downes | 7 | 18 | 12 | 48 | 21 | 9 | 1 |
| flare | 2 | 0 | 1 | 2 | 46 | 1 | 1 |
| rrlyr | 0 | 0 | 0 | 13 | 0 | 89 | 0 |
| sn | 5 | 8 | 17 | 2 | 12 | 0 | 83 |
The overall classification rate is 0.7009.
Recursive Partitioning
roz <- rpart(trform ,data=trains)
rpart.plot(roz,type=1,extra=1)
pv <- predict(roz,newdata=tests,type="class")
cm <- xtabs( ~ pv + tests$ttype)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | rrlyr | sn | |
|---|---|---|---|---|---|---|---|
| agn | 30 | 2 | 2 | 8 | 1 | 0 | 3 |
| blazar | 0 | 12 | 4 | 15 | 0 | 4 | 1 |
| cv | 2 | 6 | 88 | 37 | 5 | 3 | 6 |
| downes | 6 | 7 | 9 | 52 | 6 | 12 | 6 |
| flare | 1 | 0 | 0 | 1 | 6 | 0 | 3 |
| rrlyr | 0 | 3 | 0 | 10 | 0 | 84 | 1 |
| sn | 5 | 9 | 42 | 12 | 6 | 0 | 172 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | rrlyr | sn | |
|---|---|---|---|---|---|---|---|
| agn | 68 | 5 | 1 | 6 | 4 | 0 | 2 |
| blazar | 0 | 31 | 3 | 11 | 0 | 4 | 1 |
| cv | 5 | 15 | 61 | 27 | 21 | 3 | 3 |
| downes | 14 | 18 | 6 | 39 | 25 | 12 | 3 |
| flare | 2 | 0 | 0 | 1 | 25 | 0 | 2 |
| rrlyr | 0 | 8 | 0 | 7 | 0 | 82 | 1 |
| sn | 11 | 23 | 29 | 9 | 25 | 0 | 90 |
The overall classification rate is 0.651.
Support Vector Machines
Use the default choice of setting from the kernlab R package for this:
svmod <- ksvm(trform, data=trains)Using automatic sigma estimation (sigest) for RBF or laplace kernel pv <- predict(svmod, tests)
cm <- xtabs( ~ pv + tests$ttype)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | rrlyr | sn | |
|---|---|---|---|---|---|---|---|
| agn | 37 | 2 | 0 | 11 | 0 | 0 | 6 |
| blazar | 0 | 21 | 3 | 7 | 0 | 0 | 0 |
| cv | 2 | 8 | 108 | 36 | 2 | 2 | 17 |
| downes | 1 | 2 | 9 | 69 | 4 | 8 | 1 |
| flare | 0 | 0 | 0 | 1 | 12 | 0 | 1 |
| rrlyr | 0 | 2 | 0 | 7 | 0 | 93 | 0 |
| sn | 4 | 4 | 25 | 4 | 6 | 0 | 167 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | rrlyr | sn | |
|---|---|---|---|---|---|---|---|
| agn | 84 | 5 | 0 | 8 | 0 | 0 | 3 |
| blazar | 0 | 54 | 2 | 5 | 0 | 0 | 0 |
| cv | 5 | 21 | 74 | 27 | 8 | 2 | 9 |
| downes | 2 | 5 | 6 | 51 | 17 | 8 | 1 |
| flare | 0 | 0 | 0 | 1 | 50 | 0 | 1 |
| rrlyr | 0 | 5 | 0 | 5 | 0 | 90 | 0 |
| sn | 9 | 10 | 17 | 3 | 25 | 0 | 87 |
The overall classification rate is 0.7434.
Neural Net
Use the multinom() function from the nnet R package. Might work better with some scaling.
svmod <- multinom(trform, data=trains, trace=FALSE, maxit=1000, decay=5e-4)
pv <- predict(svmod, tests)
cm <- xtabs( ~ pv + tests$ttype)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | rrlyr | sn | |
|---|---|---|---|---|---|---|---|
| agn | 35 | 1 | 1 | 5 | 0 | 0 | 3 |
| blazar | 0 | 22 | 4 | 11 | 0 | 0 | 3 |
| cv | 3 | 5 | 93 | 26 | 0 | 1 | 18 |
| downes | 3 | 8 | 23 | 80 | 6 | 8 | 5 |
| flare | 0 | 0 | 0 | 2 | 11 | 1 | 4 |
| rrlyr | 0 | 0 | 0 | 7 | 2 | 93 | 0 |
| sn | 3 | 3 | 24 | 4 | 5 | 0 | 159 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | rrlyr | sn | |
|---|---|---|---|---|---|---|---|
| agn | 80 | 3 | 1 | 4 | 0 | 0 | 2 |
| blazar | 0 | 56 | 3 | 8 | 0 | 0 | 2 |
| cv | 7 | 13 | 64 | 19 | 0 | 1 | 9 |
| downes | 7 | 21 | 16 | 59 | 25 | 8 | 3 |
| flare | 0 | 0 | 0 | 1 | 46 | 1 | 2 |
| rrlyr | 0 | 0 | 0 | 5 | 8 | 90 | 0 |
| sn | 7 | 8 | 17 | 3 | 21 | 0 | 83 |
The overall classification rate is 0.7229.
Random Forest
Use the randomForest package with the default settings:
tallform <- as.formula(paste("ttype ~",tpredform))
fmod <- randomForest(tallform, data=na.omit(trains))
pv <- predict(fmod, newdata=na.omit(tests))
cm <- xtabs( ~ pv + na.omit(tests)$ttype)This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | rrlyr | sn | |
|---|---|---|---|---|---|---|---|
| agn | 35 | 2 | 1 | 5 | 0 | 0 | 3 |
| blazar | 0 | 22 | 3 | 5 | 0 | 0 | 0 |
| cv | 3 | 7 | 94 | 28 | 2 | 1 | 16 |
| downes | 2 | 3 | 18 | 78 | 4 | 5 | 5 |
| flare | 0 | 0 | 0 | 5 | 13 | 1 | 0 |
| rrlyr | 0 | 2 | 0 | 6 | 0 | 96 | 0 |
| sn | 4 | 3 | 29 | 8 | 5 | 0 | 168 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | rrlyr | sn | |
|---|---|---|---|---|---|---|---|
| agn | 80 | 5 | 1 | 4 | 0 | 0 | 2 |
| blazar | 0 | 56 | 2 | 4 | 0 | 0 | 0 |
| cv | 7 | 18 | 65 | 21 | 8 | 1 | 8 |
| downes | 5 | 8 | 12 | 58 | 17 | 5 | 3 |
| flare | 0 | 0 | 0 | 4 | 54 | 1 | 0 |
| rrlyr | 0 | 5 | 0 | 4 | 0 | 93 | 0 |
| sn | 9 | 8 | 20 | 6 | 21 | 0 | 88 |
The overall classification rate is 0.7419.
Heirarchical Classification
First we classify into transient and non-variable. The cases which are classified as transient are then classified into type of transient. The transient classification here is different from the one above in the data used. Above, all the data are known to be transients whereas here some cases from the non-variable set will have been classified as transient at the first stage.
LDA
ldamod <- lda(vtform ,data=trains)
pv <- predict(ldamod, tests)
utests <- subset(tests, pv$class != 'nv')
pvt <- predict(ldamod, trains)
utrains <- subset(trains, pvt$class != 'nv')
ldamod <- lda(trform, data=utrains)
predc <- as.character(pv$class)
predc[predc != 'nv'] <- as.character(predict(ldamod, utests)$class)
cm <- xtabs( ~ predc + as.character(tests$type))This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 28 | 2 | 3 | 2 | 2 | 2 | 0 | 8 |
| blazar | 0 | 24 | 6 | 12 | 0 | 0 | 0 | 5 |
| cv | 3 | 4 | 84 | 25 | 0 | 7 | 0 | 17 |
| downes | 0 | 6 | 22 | 46 | 0 | 7 | 11 | 1 |
| flare | 0 | 0 | 1 | 2 | 3 | 8 | 1 | 0 |
| nv | 12 | 0 | 9 | 30 | 16 | 523 | 0 | 17 |
| rrlyr | 0 | 0 | 0 | 17 | 0 | 2 | 91 | 0 |
| sn | 1 | 3 | 20 | 1 | 3 | 9 | 0 | 144 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 64 | 5 | 2 | 1 | 8 | 0 | 0 | 4 |
| blazar | 0 | 62 | 4 | 9 | 0 | 0 | 0 | 3 |
| cv | 7 | 10 | 58 | 19 | 0 | 1 | 0 | 9 |
| downes | 0 | 15 | 15 | 34 | 0 | 1 | 11 | 1 |
| flare | 0 | 0 | 1 | 1 | 12 | 1 | 1 | 0 |
| nv | 27 | 0 | 6 | 22 | 67 | 94 | 0 | 9 |
| rrlyr | 0 | 0 | 0 | 13 | 0 | 0 | 88 | 0 |
| sn | 2 | 8 | 14 | 1 | 12 | 2 | 0 | 75 |
The overall classification rate is 0.7605.
RPART
roz <- rpart(vtform ,data=trains)
pv <- predict(roz, tests, type="class")
utests <- subset(tests, pv != 'nv')
pvt <- predict(roz, trains, type="class")
utrains <- subset(trains, pvt != 'nv')
roz <- rpart(trform, data=utrains)
predc <- as.character(pv)
predc[predc != 'nv'] <- as.character(predict(roz, utests,
type="class"))
predc <- factor(predc, levels=sort(levels(trains$type)))
cm <- xtabs( ~ predc + as.character(tests$type))This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 33 | 2 | 0 | 5 | 3 | 17 | 0 | 11 |
| blazar | 0 | 15 | 5 | 10 | 0 | 0 | 0 | 2 |
| cv | 1 | 6 | 83 | 35 | 0 | 2 | 1 | 5 |
| downes | 0 | 9 | 11 | 34 | 0 | 7 | 17 | 9 |
| flare | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| nv | 7 | 1 | 16 | 36 | 14 | 513 | 2 | 23 |
| rrlyr | 0 | 3 | 0 | 7 | 0 | 1 | 83 | 1 |
| sn | 3 | 3 | 30 | 8 | 7 | 18 | 0 | 141 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 75 | 5 | 0 | 4 | 12 | 3 | 0 | 6 |
| blazar | 0 | 38 | 3 | 7 | 0 | 0 | 0 | 1 |
| cv | 2 | 15 | 57 | 26 | 0 | 0 | 1 | 3 |
| downes | 0 | 23 | 8 | 25 | 0 | 1 | 17 | 5 |
| flare | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| nv | 16 | 3 | 11 | 27 | 58 | 92 | 2 | 12 |
| rrlyr | 0 | 8 | 0 | 5 | 0 | 0 | 81 | 1 |
| sn | 7 | 8 | 21 | 6 | 29 | 3 | 0 | 73 |
The overall classification rate is 0.7274.
SVM
svmod <- ksvm(vtform ,data=trains)Using automatic sigma estimation (sigest) for RBF or laplace kernel pv <- predict(svmod, tests)
utests <- subset(tests, pv != 'nv')
pvt <- predict(svmod, trains)
utrains <- subset(trains, pvt != 'nv')
svmod <- ksvm(trform, data=utrains)Using automatic sigma estimation (sigest) for RBF or laplace kernel predc <- as.character(pv)
predc[predc != 'nv'] <- as.character(predict(svmod, utests))
predc <- factor(predc, levels=sort(levels(trains$type)))
cm <- xtabs( ~ predc + as.character(tests$type))This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 32 | 2 | 0 | 6 | 0 | 5 | 0 | 3 |
| blazar | 0 | 23 | 3 | 8 | 0 | 0 | 0 | 0 |
| cv | 3 | 6 | 103 | 35 | 2 | 6 | 2 | 13 |
| downes | 0 | 1 | 9 | 47 | 0 | 7 | 8 | 1 |
| flare | 0 | 0 | 0 | 1 | 8 | 4 | 0 | 1 |
| nv | 7 | 0 | 6 | 28 | 11 | 524 | 0 | 13 |
| rrlyr | 0 | 2 | 0 | 6 | 0 | 1 | 93 | 0 |
| sn | 2 | 5 | 24 | 4 | 3 | 11 | 0 | 161 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 73 | 5 | 0 | 4 | 0 | 1 | 0 | 2 |
| blazar | 0 | 59 | 2 | 6 | 0 | 0 | 0 | 0 |
| cv | 7 | 15 | 71 | 26 | 8 | 1 | 2 | 7 |
| downes | 0 | 3 | 6 | 35 | 0 | 1 | 8 | 1 |
| flare | 0 | 0 | 0 | 1 | 33 | 1 | 0 | 1 |
| nv | 16 | 0 | 4 | 21 | 46 | 94 | 0 | 7 |
| rrlyr | 0 | 5 | 0 | 4 | 0 | 0 | 90 | 0 |
| sn | 5 | 13 | 17 | 3 | 12 | 2 | 0 | 84 |
The overall classification rate is 0.7992.
NNET
svmod <- multinom(vtform, data=trains, trace=FALSE, maxit=1000, decay=5e-4)
pv <- predict(svmod, tests)
utests <- subset(tests, pv != 'nv')
pvt <- predict(svmod, trains)
utrains <- subset(trains, pvt != 'nv')
svmod <- multinom(trform, data=trains, trace=FALSE, maxit=1000, decay=5e-4)
predc <- as.character(pv)
predc[predc != 'nv'] <- as.character(predict(svmod, utests))
predc <- factor(predc, levels=sort(levels(trains$type)))
cm <- xtabs( ~ predc + as.character(tests$type))This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 31 | 1 | 0 | 3 | 0 | 2 | 0 | 3 |
| blazar | 0 | 22 | 4 | 11 | 0 | 1 | 0 | 3 |
| cv | 3 | 5 | 91 | 24 | 0 | 6 | 1 | 18 |
| downes | 2 | 8 | 23 | 64 | 0 | 13 | 8 | 5 |
| flare | 0 | 0 | 0 | 0 | 6 | 6 | 1 | 3 |
| nv | 7 | 0 | 6 | 24 | 15 | 520 | 0 | 14 |
| rrlyr | 0 | 0 | 0 | 6 | 0 | 0 | 93 | 0 |
| sn | 1 | 3 | 21 | 3 | 3 | 10 | 0 | 146 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 70 | 3 | 0 | 2 | 0 | 0 | 0 | 2 |
| blazar | 0 | 56 | 3 | 8 | 0 | 0 | 0 | 2 |
| cv | 7 | 13 | 63 | 18 | 0 | 1 | 1 | 9 |
| downes | 5 | 21 | 16 | 47 | 0 | 2 | 8 | 3 |
| flare | 0 | 0 | 0 | 0 | 25 | 1 | 1 | 2 |
| nv | 16 | 0 | 4 | 18 | 62 | 93 | 0 | 7 |
| rrlyr | 0 | 0 | 0 | 4 | 0 | 0 | 90 | 0 |
| sn | 2 | 8 | 14 | 2 | 12 | 2 | 0 | 76 |
The overall classification rate is 0.7847.
Random Forest
tallform <- as.formula(paste("vtype ~",tpredform))
svmod <- randomForest(tallform, data=trains)
pv <- predict(svmod, tests)
utests <- subset(tests, pv != 'nv')
pvt <- predict(svmod, trains)
utrains <- subset(trains, pvt != 'nv')
tallform <- as.formula(paste("ttype ~",tpredform))
svmod <- randomForest(tallform, data=na.omit(trains))
predc <- as.character(pv)
predc[predc != 'nv'] <- as.character(predict(svmod, utests))
predc <- factor(predc, levels=sort(levels(trains$type)))
cm <- xtabs( ~ predc + as.character(tests$type))This table shows the predicted type in the rows by the actual type in the columns.
print(xtable(cm,digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 31 | 2 | 1 | 4 | 0 | 4 | 0 | 2 |
| blazar | 0 | 22 | 3 | 7 | 0 | 1 | 0 | 0 |
| cv | 3 | 8 | 92 | 28 | 0 | 2 | 1 | 12 |
| downes | 0 | 2 | 17 | 57 | 0 | 11 | 8 | 3 |
| flare | 0 | 0 | 0 | 2 | 10 | 5 | 1 | 0 |
| nv | 7 | 0 | 5 | 23 | 11 | 520 | 0 | 11 |
| rrlyr | 0 | 2 | 0 | 6 | 0 | 1 | 93 | 0 |
| sn | 3 | 3 | 27 | 8 | 3 | 14 | 0 | 164 |
Same as above but now expressed a percentage within each column:
print(xtable(round(100*prop.table(cm, 2)),digits=0,caption="Actual"),type="html",caption.placement="top")| agn | blazar | cv | downes | flare | nv | rrlyr | sn | |
|---|---|---|---|---|---|---|---|---|
| agn | 70 | 5 | 1 | 3 | 0 | 1 | 0 | 1 |
| blazar | 0 | 56 | 2 | 5 | 0 | 0 | 0 | 0 |
| cv | 7 | 21 | 63 | 21 | 0 | 0 | 1 | 6 |
| downes | 0 | 5 | 12 | 42 | 0 | 2 | 8 | 2 |
| flare | 0 | 0 | 0 | 1 | 42 | 1 | 1 | 0 |
| nv | 16 | 0 | 3 | 17 | 46 | 93 | 0 | 6 |
| rrlyr | 0 | 5 | 0 | 4 | 0 | 0 | 90 | 0 |
| sn | 7 | 8 | 19 | 6 | 12 | 3 | 0 | 85 |
The overall classification rate is 0.7976.
Summary of results
Summary of percentage classification rates across tests:
print(xtable(100*cmat,digits=2),type="html")| LDA | RPart | SVM | NN | Forest | |
|---|---|---|---|---|---|
| All | 75.97 | 71.94 | 80.24 | 79.60 | 80.48 |
| TranNoTran | 90.40 | 88.39 | 92.02 | 91.61 | 91.77 |
| Tranonly | 70.09 | 65.10 | 74.34 | 72.29 | 74.19 |
| Heirarch | 76.05 | 72.74 | 79.92 | 78.47 | 79.76 |
cmaternGPR <- cmat