inst/doc/f_Multivariate_Approach_Prediction_Model.R

In mand: Multivariate Analysis for Neuroimaging Data

## ----echo=FALSE, error=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
options(warn=-1)
knitr::opts_chunk$set(eval = FALSE)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  if(!require("mand")) install.packages("mand")

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  library(mand)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  fit113 = msma(SB1, Z=Z, comp=2, lambdaX=0.075, muX=0.5)
#  Ss = fit113$ssX
#  colnames(Ss) = paste("c", c(1:ncol(Ss)), sep="")
#  swdata113 = data.frame(
#  Z = as.factor(ifelse(Z == 1, "Y", "N")), Ss)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  glmfit = glm(Z~., data=swdata113, family=binomial)
#  summary(glmfit)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  test = predict(glmfit, type="response")>=0.5

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (err.table = table(swdata113$Z, test))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  1 - sum(diag(err.table)) / sum(err.table)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  t(apply(err.table, 1, function(x) x / sum(x)))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  matrix(
#  c("specificity", "false positive rate",
#  "false negative rate","sensitivity")
#  , ncol=2)

## ----fig.width = 5, fig.height = 3.5----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  x = seq(min(swdata113$c1), max(swdata113$c1), length = 30)
#  y = seq(min(swdata113$c2), max(swdata113$c2),
#  length = length(x))
#  prob = function(x, y) 1/(1+exp(-predict(glmfit,
#  newdata=data.frame(c1=x, c2=y))))
#  z = outer(x, y, prob)
#  filled.contour(x,y,z, xlab="Component 1", ylab="Component 2")

## ----message=FALSE----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  library(e1071)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  set.seed(1)
#  tuneSVM = tune(svm, Z~., data=swdata113,
#  ranges = list(gamma = 2^(0:2), cost = c(4, 6, 8)),
#  tunecontrol = tune.control(cross = nrow(swdata113)))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  summary(tuneSVM)

## ----fig.width = 5, fig.height = 3.5----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  plot(tuneSVM, color.palette = heat.colors)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  bestGamma = tuneSVM$best.parameters$gamma
#  bestC = tuneSVM$best.parameters$cost

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  set.seed(1)
#  svmfit = svm(Z~., data=swdata113,
#  cost = bestC, gamma = bestGamma,
#  probability=TRUE, kernel="radial", cross=nrow(swdata113))
#  summary(svmfit)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  pred = predict(svmfit, newdata=swdata113, probability=TRUE,
#  decision.values=TRUE)
#  (err.table = table(swdata113$Z, pred))
#  1 - sum(diag(err.table)) / sum(err.table)

## ----fig.width = 5, fig.height = 3.5----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  plot(svmfit, swdata113, c2~c1)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  opt11 = optparasearch(SB1, Z=Z, comp=20,
#  search.method = "regparaonly", criterion="BIC")
#  (fit311 = msma(SB1, Z=Z,
#  comp=opt11$optncomp, lambdaX=opt11$optlambdaX))
#  Ss = fit311$ssX
#  colnames(Ss) = paste("c", c(1:ncol(Ss)), sep="")
#  swdata311 = data.frame(
#  Z = as.factor(ifelse(Z == 1, "Y", "N")), Ss)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  if(!require("rpart.plot")) install.packages("rpart.plot")
#  library(rpart)
#  library(rpart.plot)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  set.seed(1)
#  (treefit = rpart(Z~., data=swdata311,
#  control = rpart.control(minsplit = 4)))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  prp(treefit, type=4, extra=1, faclen=0, nn=TRUE)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  printcp(treefit)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  plotcp(treefit)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (treefit1 = prune(treefit, cp=0.05))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (treefit2 = snip.rpart(treefit, 7))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  prp(treefit2, type=4, extra=1, faclen=0, nn=TRUE)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  pred = predict(treefit2, type="class")

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (err.table = table(swdata311$Z, pred))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  1 - sum(diag(err.table))/sum(err.table)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  swdata3112 = head(swdata311)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  set.seed(1)
#  (idrand = sample(1:6, replace=TRUE))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (bsample = swdata3112[idrand, 1:4])

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (oobsample = swdata3112[!(1:nrow(swdata3112) %in%
#  unique(idrand)), 1:4])

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  library(randomForest)
#  library(e1071)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  set.seed(1)
#  tuneRF = tune(randomForest, Z~., data=swdata311,
#  ranges = list(mtry = c(4,6,8), ntree = c(300, 500, 1000),
#  nodesize= c(1,2,3)),
#  tunecontrol = tune.control(cross = nrow(swdata311)))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  summary(tuneRF)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  bestmtry = tuneRF$best.parameters$mtry
#  bestntree = tuneRF$best.parameters$ntree
#  bestnodesize = tuneRF$best.parameters$nodesize

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  set.seed(1)
#  (rffit = randomForest(Z~., data=swdata311, proximity=TRUE,
#  mtry = bestmtry, ntree = bestntree, nodesize=bestnodesize))

## ----fig.width = 5, fig.height = 5------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  varImpPlot(rffit)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  Q = fit311$wbX[[1]][,6]
#  outstat1 = rec(Q, img1$imagedim, B=B1, mask=img1$brainpos)
#  outstat2 = -outstat1
#  coat(template, outstat2)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  atlastable(tmpatlas, outstat2, atlasdataset)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  pred = predict(rffit, type="class")

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (err.table = table(swdata311$Z, pred))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  1 - sum(diag(err.table))/sum(err.table)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  t(apply(err.table, 1, function(x) x / sum(x)))

## ----fig.width = 6, fig.height = 3.5----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  par(mfrow=c(1,2))
#  partialPlot(rffit, swdata311, c1)
#  partialPlot(rffit, swdata311, c2)

## ----fig.width = 4.5, fig.height = 3----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  par(mfrow=c(1,1), mar=c(3,3,3,8))
#  z = rffit$proximity
#  n = nrow(swdata311)
#  filled.contour(x=1:n, y=1:n, z=z, color = terrain.colors)

## ----fig.width = 4.5, fig.height = 3----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  par(mfrow=c(1,1), mar=c(4,3,2,2))
#  MDSplot(rffit, factor(swdata311$Z),
#  pch=as.numeric(swdata311$Z)-1)

## ----fig.width = 3.5, fig.height = 3----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  par(mfrow=c(1,1), mar=c(4,3,2,2))
#  plot(randomForest::outlier(rffit), type="h",
#  col= as.numeric(swdata311$Z))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  img2 = simbrain(baseimg = baseimg, diffimg = diffimg2,
#  sdevimg=sdevimg, mask=mask, n0=500, c1=0.01, sd1=0.1,
#  zeromask=FALSE, seed=2)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  testZ = as.factor(ifelse(img2$Z == 1, "Y", "N"))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  SB2 = basisprod(img2$S, B1)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  ptest113 = ptest(object=fit113, Z=Z, newdata=SB2,
#  testZ=testZ, regmethod = "glm")

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  summary(ptest113$trainout$finalModel)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  ptest113$trainout

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  ptest113$predcnfmat

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  ptest112 = ptest(object=fit112, Z=Z, newdata=SB2, testZ=testZ,
#  regmethod = "glm")
#  ptest112$predcnfmat$overall["Accuracy"]

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  ptest311 = ptest(object=img1$S, Z=Z, newdata=img2$S,
#  testZ=testZ, regmethod = "glm")
#  ptest312 = ptest(object=SB1, Z=Z, newdata=SB2, testZ=testZ,
#  regmethod = "glm")

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  ptest311$predcnfmat$overall["Accuracy"]
#  ptest312$predcnfmat$overall["Accuracy"]

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  mus = seq(0, 1, by=0.25)
#  comps = c(1,2,5,10)
#  
#  paramtest1 = lapply(comps, function(c1){lapply(mus,
#  function(mu1){
#  tmpfit = msma(SB1, Z=Z, comp=c1, lambdaX=0.075, muX=mu1)
#  tmpptest = ptest(object=tmpfit, Z=Z, newdata=SB2,
#  testZ=testZ, regmethod = "glm")
#  tmpptest$predcnfmat
#  })})
#  
#  out1 = do.call(rbind, lapply(paramtest1, function(x)
#  do.call(cbind, lapply(x,
#  function(y)y$overall["Accuracy"]))))
#  rownames(out1)=comps; colnames(out1)=mus
#  

## ----results='asis'---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  kable(out1, "latex", booktabs = T)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  svmgrid = expand.grid(sigma = c(0.001, 0.025, 0.05),
#  C = c(0.5, 0.75, 1))
#  ptest211 = ptest(object=fit311, Z=Z, newdata=SB2, testZ=testZ,
#  regmethod = "svmRadial", metric="ROC",
#  param=svmgrid)
#  ptest211$trainout

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  ptest211$predcnfmat

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  treegrid = NULL
#  ptest212 = ptest(object=fit311, Z=Z, newdata=SB2, testZ=testZ,
#  regmethod = "rpart", metric="ROC",
#  param=treegrid)
#  ptest212$trainout

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  ptest212$predcnfmat

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  rfgrid  =  NULL
#  ptest213 = ptest(object=fit311, Z=Z, newdata=SB2, testZ=testZ,
#  regmethod = "rf", metric="ROC",
#  param=rfgrid)
#  ptest213$trainout

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  ptest213$predcnfmat

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  layers0 = c(1, 5, 10); layers1 = c(0, 1, 5, 10)
#  rate0 = c(0, 0.25, 0.5, 0.75)
#  activation=c("relu", "sigmoid", "tanh", "softrelu")
#  mxnet.params = expand.grid(layer1=layers0, layer2=layers1,
#  layer3=0, learning.rate=0.1, momentum=0.9, dropout=0,
#  activation=activation[3])

## ----eval=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  ptest215 = ptest(object=fit113, Z=Z, newdata=SB2, testZ=testZ,
#  regmethod = "mxnet", metric="Accuracy",
#  param=mxnet.params)
#  ptest215$trainout