20211506_Data

options(repos = c(CRAN = "https://cran.r-project.org/"))
install.packages("GGally")

## 'C:/Users/billd/AppData/Local/R/win-library/4.3'의 위치에 패키지(들)을 설치합니다.
## (왜냐하면 'lib'가 지정되지 않았기 때문입니다)

## 패키지 'GGally'를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다
## 
## 다운로드된 바이너리 패키지들은 다음의 위치에 있습니다
##  C:\Users\billd\AppData\Local\Temp\Rtmpwv0Bu8\downloaded_packages

install.packages("caret")

## 'C:/Users/billd/AppData/Local/R/win-library/4.3'의 위치에 패키지(들)을 설치합니다.
## (왜냐하면 'lib'가 지정되지 않았기 때문입니다)

## 패키지 'caret'를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

## Warning: 패키지 'caret'의 이전설치를 삭제할 수 없습니다

## Warning in file.copy(savedcopy, lib, recursive = TRUE):
## C:\Users\billd\AppData\Local\R\win-library\4.3\00LOCK\caret\libs\x64\caret.dll를
## C:\Users\billd\AppData\Local\R\win-library\4.3\caret\libs\x64\caret.dll로
## 복사하는데 문제가 발생했습니다: Permission denied

## Warning: 'caret'를 복구하였습니다

## 
## 다운로드된 바이너리 패키지들은 다음의 위치에 있습니다
##  C:\Users\billd\AppData\Local\Temp\Rtmpwv0Bu8\downloaded_packages

install.packages("partykit")

## 'C:/Users/billd/AppData/Local/R/win-library/4.3'의 위치에 패키지(들)을 설치합니다.
## (왜냐하면 'lib'가 지정되지 않았기 때문입니다)

## 패키지 'partykit'를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

## Warning: 패키지 'partykit'의 이전설치를 삭제할 수 없습니다

## Warning in file.copy(savedcopy, lib, recursive = TRUE):
## C:\Users\billd\AppData\Local\R\win-library\4.3\00LOCK\partykit\libs\x64\partykit.dll를
## C:\Users\billd\AppData\Local\R\win-library\4.3\partykit\libs\x64\partykit.dll로
## 복사하는데 문제가 발생했습니다: Permission denied

## Warning: 'partykit'를 복구하였습니다

## 
## 다운로드된 바이너리 패키지들은 다음의 위치에 있습니다
##  C:\Users\billd\AppData\Local\Temp\Rtmpwv0Bu8\downloaded_packages

install.packages("rpart")

## 'C:/Users/billd/AppData/Local/R/win-library/4.3'의 위치에 패키지(들)을 설치합니다.
## (왜냐하면 'lib'가 지정되지 않았기 때문입니다)

## 패키지 'rpart'를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다
## 
## 다운로드된 바이너리 패키지들은 다음의 위치에 있습니다
##  C:\Users\billd\AppData\Local\Temp\Rtmpwv0Bu8\downloaded_packages

install.packages("ggplot2")

## 'C:/Users/billd/AppData/Local/R/win-library/4.3'의 위치에 패키지(들)을 설치합니다.
## (왜냐하면 'lib'가 지정되지 않았기 때문입니다)

## 패키지 'ggplot2'를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다
## 
## 다운로드된 바이너리 패키지들은 다음의 위치에 있습니다
##  C:\Users\billd\AppData\Local\Temp\Rtmpwv0Bu8\downloaded_packages

install.packages("patchwork")

## 'C:/Users/billd/AppData/Local/R/win-library/4.3'의 위치에 패키지(들)을 설치합니다.
## (왜냐하면 'lib'가 지정되지 않았기 때문입니다)

## 패키지 'patchwork'를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다
## 
## 다운로드된 바이너리 패키지들은 다음의 위치에 있습니다
##  C:\Users\billd\AppData\Local\Temp\Rtmpwv0Bu8\downloaded_packages

install.packages("tidyverse")

## 'C:/Users/billd/AppData/Local/R/win-library/4.3'의 위치에 패키지(들)을 설치합니다.
## (왜냐하면 'lib'가 지정되지 않았기 때문입니다)

## 패키지 'tidyverse'를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다
## 
## 다운로드된 바이너리 패키지들은 다음의 위치에 있습니다
##  C:\Users\billd\AppData\Local\Temp\Rtmpwv0Bu8\downloaded_packages

install.packages("nnet")

## 'C:/Users/billd/AppData/Local/R/win-library/4.3'의 위치에 패키지(들)을 설치합니다.
## (왜냐하면 'lib'가 지정되지 않았기 때문입니다)

## 패키지 'nnet'를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다
## 
## 다운로드된 바이너리 패키지들은 다음의 위치에 있습니다
##  C:\Users\billd\AppData\Local\Temp\Rtmpwv0Bu8\downloaded_packages

library(GGally)

## Warning: 패키지 'GGally'는 R 버전 4.3.3에서 작성되었습니다

## 필요한 패키지를 로딩중입니다: ggplot2

## Warning: 패키지 'ggplot2'는 R 버전 4.3.3에서 작성되었습니다

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

library(caret)

## Warning: 패키지 'caret'는 R 버전 4.3.3에서 작성되었습니다

## 필요한 패키지를 로딩중입니다: lattice

library(partykit)

## Warning: 패키지 'partykit'는 R 버전 4.3.3에서 작성되었습니다

## 필요한 패키지를 로딩중입니다: grid

## 필요한 패키지를 로딩중입니다: libcoin

## Warning: 패키지 'libcoin'는 R 버전 4.3.3에서 작성되었습니다

## 필요한 패키지를 로딩중입니다: mvtnorm

## Warning: 패키지 'mvtnorm'는 R 버전 4.3.3에서 작성되었습니다

library(rpart)

## Warning: 패키지 'rpart'는 R 버전 4.3.3에서 작성되었습니다

library(ggplot2)
library(patchwork)

## Warning: 패키지 'patchwork'는 R 버전 4.3.3에서 작성되었습니다

library(tidyverse)

## Warning: 패키지 'tidyverse'는 R 버전 4.3.3에서 작성되었습니다

## Warning: 패키지 'tibble'는 R 버전 4.3.3에서 작성되었습니다

## Warning: 패키지 'tidyr'는 R 버전 4.3.3에서 작성되었습니다

## Warning: 패키지 'readr'는 R 버전 4.3.3에서 작성되었습니다

## Warning: 패키지 'purrr'는 R 버전 4.3.3에서 작성되었습니다

## Warning: 패키지 'dplyr'는 R 버전 4.3.3에서 작성되었습니다

## Warning: 패키지 'stringr'는 R 버전 4.3.3에서 작성되었습니다

## Warning: 패키지 'forcats'는 R 버전 4.3.3에서 작성되었습니다

## Warning: 패키지 'lubridate'는 R 버전 4.3.3에서 작성되었습니다

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ purrr::lift()   masks caret::lift()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(nnet)

## Warning: 패키지 'nnet'는 R 버전 4.3.3에서 작성되었습니다

install.packages("conflicted")

## 'C:/Users/billd/AppData/Local/R/win-library/4.3'의 위치에 패키지(들)을 설치합니다.
## (왜냐하면 'lib'가 지정되지 않았기 때문입니다)

## 패키지 'conflicted'를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다
## 
## 다운로드된 바이너리 패키지들은 다음의 위치에 있습니다
##  C:\Users\billd\AppData\Local\Temp\Rtmpwv0Bu8\downloaded_packages

library(conflicted)

## Warning: 패키지 'conflicted'는 R 버전 4.3.3에서 작성되었습니다

conflict_prefer("filter", "dplyr")

## [conflicted] Will prefer dplyr::filter over any other package.

conflict_prefer("lag", "dplyr")

## [conflicted] Will prefer dplyr::lag over any other package.

conflict_prefer("lift", "purrr")

## [conflicted] Will prefer purrr::lift over any other package.

set_olive <- function(inc=c(1,2)){
  require(tidyverse, quietly = TRUE)
  
  if(!("cepp" %in% .packages(all = TRUE))) install.packages("cepp") 
  require(cepp, quietly = TRUE)
  data(olive)
  
  olivex <- olive
  names(olivex) <- c("rg","ar","pt","po","st","ol","lo","ln","ac","es")
  olivex
  rg <- factor(olivex$rg)
  levels(rg) <- c("S","R","N")
  olivex$rg <- factor(rg, levels=c("S","R","N"))
  
  ar <- factor(olivex$ar)
  arf <- factor(ar, levels=levels(ar)[c(6,1,7,5,2,4,8,3,9)])
  
  x <- levels(arf) %>% 
    tolower() %>% 
    str_replace("east","e") %>% 
    str_replace("west","w") %>%
    str_replace("south","s") %>%
    str_replace("north","n") %>%
    str_replace("coast","c") %>%
    str_replace("inland","i") 
  
  y <- str_split(x,"-",simplify=TRUE)
  y1 <- str_sub(y[,1],1,4)
  y2 <- str_sub(y[,2],1,4)
  
  short <- paste0(c("S","S","S","S","R","R","N","N","N"),
                  "_", y1,"_", y2) %>% str_replace_all("_$","")
  
  levels(arf) <- short
  olivex$ar <- arf
  
  indx <- c(unique(c(inc, 3:10)))
  olivex_subset <- olivex[, indx]
  
  return(olivex_subset) 
}

# olive1, olive2 생성
olive1 <- set_olive(1)
olive2 <- set_olive(2)

# olive1과 olive2에서 테이블 생성하기
with(olive1, table(rg))

## rg
##   S   R   N 
## 323  98 151

with(olive2, table(ar))

## ar
##   S_sici   S_cala S_s_apul S_n_apul R_c_sard R_i_sard   N_umbr N_e_ligu 
##       36       56      206       25       33       65       51       50 
## N_w_ligu 
##       50

with(set_olive(), table(rg, ar))

##    ar
## rg  S_sici S_cala S_s_apul S_n_apul R_c_sard R_i_sard N_umbr N_e_ligu N_w_ligu
##   S     36     56      206       25        0        0      0        0        0
##   R      0      0        0        0       33       65      0        0        0
##   N      0      0        0        0        0        0     51       50       50

type_fix <- function(res){
  type=class(res)
  if('multinom' %in% type) type='multinom' # nnet
  if('lda' %in% type) type='lda'     # MASS
  if('qda' %in% type) type='qda'     # MASS
  if('nnet' %in% type) type='nnet'   # nnet
  if('rpart' %in% type) type='rpart' # rpart
  if('party' %in% type) type='party' # partikit
  if('svm' %in% type) type='svm'     # e1071
  if('train' %in% type) type='train' # caret
  type
}




predict_cases<-function(res,new){
  
  type = type_fix(res)
  
  p_sp =rep( '0', nrow(new) )
  
  if(type %in% c('lda','qda')) p_sp  = 
    predict(res, newdata=new)$class  
  if(type %in% c('multinom','rpart','nnet')) p_sp  =  
    predict(res, newdata=new, type="class") 
  if(type %in% c('party','svm')) p_sp  =  
    as.vector(predict(res, newdata=new))
  if(type %in% c('train')) p_sp  =  
    predict(res, newdata=new, type="raw")
  as.character(p_sp)
}



is_a_leaf <-function(fm){
  
  out = TRUE
  if( is.call(fm) ) 
    if( length(fm)>2 )
      if( (fm[[1]]=='+')|(fm[[1]]=='-') )  out=FALSE
  out      
}  



leafs <- function(res){
  v = NULL
  fx = res
  while(!is_a_leaf(fx)){
    v = c(fx[[3]],v)
    fx = fx[[2]]
  }
  c(fx,v)
}




Draw <- function(res, border=1, axis=1, np=100){
  
  type = type_fix(res)
  
  callx  =  res$call
  df     =  callx[[3]]
  fm     =  callx[[2]]  
  
  if( type %in% c('party') ){
    df0    =  res$data
    df     =  as.name('df0')
    fm     =  res$terms       }
  
  vy     =  fm[[2]]
  terms  =  leafs(fm[[3]])
  cterms =  as.character(terms)
  pos    =  which(sapply(terms, is.name))
  
  q1     =  as.name(cterms[pos][1])
  q2     =  as.name(cterms[pos][2]) 
  
  if(axis==0){
    q1     =  as.name(cterms[pos][2])
    q2     =  as.name(cterms[pos][1])  }
  
  x1     =  eval(q1, eval(df))  
  x2     =  eval(q2, eval(df))
  
  x1q    =  seq( min(x1), max(x1), l=np)
  x2q    =  seq( min(x2), max(x2), l=np)
  
  new    =  expand.grid( x=x1q,  y=x2q )
  names(new) = as.character(c(q1,q2))
  
  cp_0   =  predict_cases(res,eval(df))
  cp_1   =  as.character(eval(vy, eval(df))) 
  p_1    =  factor(cp_1)
  dfz    =  eval(df)[cp_0!=cp_1,]
  dfz$p  =  p_1[cp_0!=cp_1]
  
  cp_sp  =  predict_cases(res,new)
  p_sp   =  factor(cp_sp, levels=levels(p_1))
  new_p  =  data.frame( new, p_sp )
  
  p<- eval( bquote( 
    ggplot(.(df)) + 
      geom_point( aes(.(q1),.(q2), col=.(vy)))+
      geom_point(data= new_p, 
                 aes(.(q1),.(q2), col=p_sp),size=0.1,alpha=0.2) +
      geom_point(data=dfz,
                 aes(.(q1),.(q2),col=p ),shape=1,size=4)+
      geom_point( aes(.(q1),.(q2), col=.(vy)), size=1.4)
    #+theme_light()
  ) )
  if(border==1){
    z_sp = factor(cp_sp)
    K = length(levels(z_sp))
    if(K>1){
      zzz<-diag(K)[z_sp,]
      for(k in 1:K){
        new_z <- new_p
        new_z$z <- zzz[,k]
        p<-p + eval( bquote( geom_contour(data=new_z,
                                          aes(.(q1),.(q2),z=z), 
                                          breaks=0.5, size=0.1, col="black") 
        ))
      }
    }
  }
  p
}

#학번이 6으로 끝나므로 2. palmitoleic ( po )과 4. oleic ( ol ) 데이터 분석
roc_LDA  <- train(ar ~ ol + po, olive2, method="lda")
roc_NNET <- train(ar ~ ol + po, olive2, method="nnet",
                  maxit=300, trace=F)
roc_TREE    <- train(ar ~ ol + po, olive2, method="rpart")
roc01a_TREE <- as.party(rpart(ar ~ ol + po, olive2, cp=0.01))
roc_SVM  <- train(ar ~ ol + po, olive2, method="svmRadial")
( pm <- t(roc_TREE$bestTune))

##            1
## cp 0.1229508

x1<- table(olive2$ar, predict(roc_LDA))
x2<- table(olive2$ar, predict(roc_NNET))
x3<- table(olive2$ar, predict(roc_TREE))
x301a<- table(olive2$ar, predict(roc01a_TREE))
x4<- table(olive2$ar, predict(roc_SVM))

sum(diag(x1))/sum(x1)

## [1] 0.7185315

sum(diag(x2))/sum(x2)

## [1] 0.7062937

sum(diag(x3))/sum(x3)

## [1] 0.5594406

sum(diag(x301a))/sum(x301a)

## [1] 0.777972

sum(diag(x4))/sum(x4)

## [1] 0.7727273

p1<- Draw(roc_LDA)+labs(title="LDA")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

p2<- Draw(roc_NNET)+labs(title="NNET")
p3<- Draw(roc01a_TREE)+labs(title="TREE")
p4<- Draw(roc_SVM)+labs(title="SVM")
(p1+p2)/(p3+p4)+plot_layout(guides="collect")

#학번이 6으로 끝나므로 1. palmitic ( pt )과 5. linolenic ( ln ) 데이터 분석
roc_LDA  <- train(ar ~ pt + ln, olive2, method="lda")
roc_NNET <- train(ar ~ pt + ln, olive2, method="nnet",
                  maxit=300, trace=F)
roc_TREE    <- train(ar ~ pt + ln, olive2, method="rpart")
roc01a_TREE <- as.party(rpart(ar ~ pt + ln, olive2, cp=0.01))
roc_SVM  <- train(ar ~ pt + ln, olive2, method="svmRadial")
( pm <- t(roc_TREE$bestTune) )

##             1
## cp 0.04644809

x1<- table(olive2$ar, predict(roc_LDA))
x2<- table(olive2$ar, predict(roc_NNET))
x3<- table(olive2$ar, predict(roc_TREE))
x301a<- table(olive2$ar, predict(roc01a_TREE))
x4<- table(olive2$ar, predict(roc_SVM))

sum(diag(x1))/sum(x1)

## [1] 0.6416084

sum(diag(x2))/sum(x2)

## [1] 0.6555944

sum(diag(x3))/sum(x3)

## [1] 0.5541958

sum(diag(x301a))/sum(x301a)

## [1] 0.736014

sum(diag(x4))/sum(x4)

## [1] 0.708042

p1<- Draw(roc_LDA)+labs(title="LDA")
p2<- Draw(roc_NNET)+labs(title="NNET")
p3<- Draw(roc01a_TREE)+labs(title="TREE")
p4<- Draw(roc_SVM)+labs(title="SVM")
(p1+p2)/(p3+p4)+plot_layout(guides="collect")

20211506_Data_Touring

김도영

2024-04-02