Random Forests regression analysis

RF 회귀분석 소개하고 내용 간단히 설명.. 

--------

# Random Forests regression analysis

# by Alan
# 24 Mar 2013

library(randomForest) # 일단 라이브러리가 설치 되어 있어야 한다. 없으면 install.packages("randomForest") 로 설치 

set.seed(1978) # 어떤 환경에서도 동일한 결과를 얻기 위해 난수 발생기에 seed 설정. 

setwd("~/Dropbox/BayesianBird/") # 작업 디렉토리 설정 

#regression by random forest

data.full <- read.csv("LKM_RFvariable5.csv") # 일단 데이터를 읽고, 

str(data.full) # 데이터 구조를 보자 

colnames(data.full) # 역시 데이터 컬럼 목록 출력 

# Factorize (very important!!!) # 이 부분이 중요. 현재 데이터에서 많은 변수가 사실 범주형 변수인데 그냥 integer 로 인식이 된다. 숫자로 코딩이 된 경우 컴퓨터가 자동으로 이게 연속 변수인지 명목 변수인지 알 수 없기 때문에 어쩔 수 없이 별도로 확인해야 함 

data.full$kizuki_1 <- factor(data.full$kizuki_1) # 여기선 Kizuki 의 P/A 자료 

data.full$major_1 <- factor(data.full$major_1) # 여기선 Major 의 P/A 자료 

data.full$Forest_management <- factor(data.full$Forest_management) 

data.full$NumofVisitor2 <- factor(data.full$NumofVisitor2)

data.full$Urbanization <- factor(data.full$Urbanization) # 등등 원래 변수 목록을 펼쳐놓고 Factor 인데 Factor 로 들어가지 않은 변수들 (대부분 숫자로 코딩된 경우)을 다시 팩터화 해 줌. 이렇게 하지 않으면 올바른 정보가 분석에 들어가지 못함. 

# Response and meta info. variables # 이후 작업을 편하게 하기 위해 컬럼 이름을 일일히 써줌. 전체 공변량을 response, plot data, environmental data 세 종류로 나눴다. 이후 environmental variable 만 가지고 맵을 만들어야 할 수도 있기 때문에 시험삼아 해 봄. 그래도 분석의 중심은 plot 수준에서 측정한 자료들이 되어야 할 것이다.  

response.vars <- c("kname", # 

                   "kizuki_1", # Factor: present(1)/absent(0)

                   "major_1", #  Factor: present(1)/absent(0)

                   "kizuki", # Float : Mean kizuki count for the 4 sampling events

                   "major") # Float: Mean major count for the 4 sampling events

# Plot-level covariates

plot.vars <- c("LogArea", # Logarithm of area size (log m2) 

               "Vegetype", # Factor: deciduous(D) / mixed (M) 

               "Habitattype", # Factor: forest (F) / ornamental trees (O)

               "Greencover1000", # Float: NDVI gradient [0,1]

               "Forest_management", # Factor: managed(1)/unmanaged(0)

               "siteNDVI", # Float: plot NDVI [0,1]

               "NumofVisitor1", # Factor: small(S) / medium(M) / large (L) 

               "NumofVisitor2", # Factor: small (0) / large (1) ?? 

               "Vege1", # Float: Normalized vegetation score1

               "Vege2"# Float: Normalized vegetation score1

               )  

# Environmental covariates

env.vars <- c("Matrix_type", # Factor (Vegetated/Agri/Urban)

              "Urbanization", # Ordered factor (1/2/3/4)

              "Road_density_500m", "Road_density_1km", "Road_density_2km", # Float: road density (%)

              "distroad_min", "distroad_avg", # Float: distance to the nearest road (m)

              "dist_f10ha", "dist_f50ha","dist_f100ha", # Float: distance to the nearest forest (m)

              "vegear_AR50", "vegear_AR40", "vegear_AR30", "vegear_AR20", "vegear_AR10", # ?? 

              "B10_AR", "B20_AR", "B30_AR", "B40_AR", "B50_AR", # Float: Buffer vegetation cover (%)

              "Greencover500", "Greencover400","Greencover300", "Greencover200", "Greencover100", # Float: Buffer NDVI [0,1]

              "dPC_PCA" # Float:  dPC = dPCintra + dPCflux + dPCconnector

              ) 

# 별도의 데이터 셋을 세 개 만들어줌. 

data.onlyEnv <- data.full[, c("kizuki", env.vars)] # 오직 환경변수만 이용하여 분석하기 위해

data.onlyPlot <- data.full[, c("kizuki", plot.vars)] # 오직 사이트 변수만 이용하여 분석하기 위해 

data.both <- data.full[, c("kizuki", env.vars, plot.vars)] # 두 종류의 데이터를 한꺼번에 이용하기 위해 

# 각 경우에 대해 별도로 RF 회귀 분석 수행. 3 가지 모형을 만들었다. 

rf.onlyEnv <- randomForest(kizuki ~ ., data= data.onlyEnv, do.trace=250, ntree=1000, mtry=4, nodesize=2, proximity=T, importance=T, corr.bias=T)

rf.onlyPlot <- randomForest(kizuki ~ ., data= data.onlyPlot, do.trace=250, ntree=1000, mtry=4, nodesize=2, proximity=T, importance=T, corr.bias=T)

rf.both <-  randomForest(kizuki ~ ., data= data.both, do.trace=250, ntree=1000, mtry=4, proximity=T, importance=T, nodesize=2, corr.bias=T)

### Tune RF    # 꼭 필요한 내용은 아니고, 최적의 mtry 파라미터 (한 트리에서 임의로 선택할 predictor 의 숫자) 를 찾기 위해. 꼭 필요한 작업은 아니다. 지정하지 않을 경우 mtry 는 predictor 갯수의 1/3 로 지정된다. 

rf.both.tune <- tuneRF(x= data.both[,-1], y=data.both$kizuki, mtryStart=3, ntreeTry=1000, trace=T, doBest=T, plot=T)

# Use mtry=12 based on the tuning # 일단 자동 탐색을 통해선 mtry = 12 로 하는 게 좋다는 결과 얻어짐 

rf.both <-  randomForest(kizuki ~ ., data= data.both, do.trace=250, ntree=3000, mtry=12, proximity=T, importance=T, nodesize=2, corr.bias=T)

#variable importance measures

# 변수 중요성 각가 플롯해 봄. 세 경우 공통적으로 log area 와 연결성 지수 (dPC: 여러 공간 범위의 dPC 를 주성분 분석 이용해 차원 축소한 것) 이 중요하게 나옴. 이전의 분석에선 선형 모형에서 중요하다고 나타나는 변수인 log area 가 중요하지 않게 나타났는데 이번엔 선형 모형과 합치되는 결과 나옴. 다만 앞서의 경우도 설명은 가능함. 사실 패치 크기 자체를 새들이 좋아하는 것은 아니고 패치 크기가 크면 새들이 좋아하는 다른 것들 (숲, 먹이, 도로에서 멀고 등) 도 같이 좋아지는 경향이 있기 때문에 (변수 간의 공선성이 존재하기 때문에) log area  가 마치 새가 많이 나오는 것을 설명하는 것처럼 보일 뿐이다. 몇 가지 지적해야 하는데, 일단 선형 모형에서 log area 와 공선성이 있기 때문에 변수 선택 과정에서 탈락하는 변수들에 아무 정보가 없느냐? 그렇지 않다. 다른 변수와 상관관계로 설명되는 부분을 뺀 나머지는 그 변수 자체의 독특한 정보고, 선형 모형에서 변수를 뺀다는 건 공선성이 없어야 한다는 가정 때문에 어쩔 수는 없지만 결국 자료에 들어있는 정보를 버리는 결과가 된다. 공선성이 존재하는 경우에도 변수를 포함시킬 수 있는 방법이 있다면 원래 가지고 있던 정보를 충분히 활용할 수 있고, RF는 그것을 수행하는 한 가지 방법. 그래서 RF 등 방법론에선 independent variable 이란 말을 쓰지 않고 covariates 나 predictor 란 표현을 쓴다. indep. var. 이란 말 자체에 다른 (독립) 변수들로부터 통계적으로 독립이어야 한다는 중요한 가정이 들어있는 것임. 

# 그래서 선형 모형에서 log area 만 남고 다 탈락 하는 것은 사실 슬픈 일이다. 우리가 가진 자료가 여러가지 있지만 공선성 문제로 인해 모든 변수들과 다 약간씩 상관관계가 있어 마치 proxvy 처럼 기능하는  log area 하나만 남게 되는 것인데, 이건 한 편으로는 size 와 diversity 간의 관계에 대한 전통적 생태학을 떠올리게 해서 reassuring 하긴 하지만 사실 새로운 것을 아무것도 얘기해 주지 않는다. 패치가 크면 새가 많아란 얘기는 지난 50년 간 해 온 얘기이기 때문에. RF 등 보다 현대적인 분석의 힘은 공선성이라는 가정으로 부터 분석을 해방시켜서 그 전에는 보이지 않던 미묘한 (고계) 상관관계를 볼 수 있게 해준다는 점..  그리고 정보를 버리지 않게 된다는 점. 결정적으로 예측력이 높다. 거의 99%의 논문에서 RF, SVM, NN, BART 등의 현대적인 방법론이 LM을 outperform 한다.,. 

varImpPlot(rf.onlyPlot)

varImpPlot(rf.onlyEnv)

varImpPlot(rf.both)

# 선형모형과 비교해 보기 위해 GLM 을 간단히 구축 

glm.both <- glm(kizuki ~ ., data=data.both) # 

# Stepwise variable selection 

# 자동 변수 선택 기능을 이용해 변수를 줄였다. 

both.glm.optimized <- step(glm.both,trace=T,  direction="both")

summary(both.glm.optimized) # 37개의 공변량 중 19개가 독립 변수로 살아남았다. 근데 요약 결과 보면 p-value 기준으로 탈락시켜야 할 변수가 아직 많이 남아 있음 (e.g. Urbanization4) 

both.glm.optimized$anova

plot(both.glm.optimized) # Diagnostic plot

##### 추가 - P-value 기준으로 2차 선택 

both.glm.optimized$R # Done already 

p.values <- summary(both.glm.optimized)$coefficients[,4]

p.values <- p.values[order(p.values, decreasing=T)]

#> p.values

 #    Urbanization2      Urbanization1      Urbanization4      Urbanization3       distroad_avg 

 #    0.9621202707       0.7999041617       0.7849329836       0.4052025878       0.2987590180 

 #   distroad_min      Greencover300            dPC_PCA             B10_AR              Vege1 

 #  0.2355416020       0.2168750390       0.1564132782       0.1458322208       0.1429523637 

 # Road_density_2km          VegetypeM Forest_management1        vegear_AR10         dist_f10ha 

 #     0.1142617807       0.0731909236       0.0629673645       0.0616547843       0.0579060275 

 # Road_density_1km        vegear_AR40        vegear_AR30      Greencover200              Vege2 

 #     0.0556675908       0.0382576639       0.0312023774       0.0291330865       0.0266012788 

 #           B30_AR             B20_AR      Greencover400  Road_density_500m      Greencover500 

 #     0.0247461292       0.0169531387       0.0024097807       0.0022316875       0.0021394711 

 #    Greencover100 

  #    0.0005496224 

barplot(names.arg=names(p.values),height= p.values, ylim=c(0, 1), cex.names=0.5) # P-values 

variables.selected.2 <- labels(p.values[which(p.values < 0.05)]) # 16 variables don't pass the cut-off value (a=0.16)

variables.selected.2[variables.selected.2=="VegetypeM"] <- "Vegetype"

variables.selected.2[variables.selected.2=="Forest_management1"] <- "Forest_management"

# 10개 가 a=0.05 하에서 살아남음 

# [1] "vegear_AR40"       "vegear_AR30"       "Greencover200"     "Vege2"            

# [5] "B30_AR"            "B20_AR"            "Greencover400"     "Road_density_500m"

# [9] "Greencover500"     "Greencover100"  

both.glm.optimized.2 <- glm(kizuki ~ . -1, data=subset(data.both, select=c("kizuki", variables.selected.2)))

both.glm.optimized.2

summary(both.glm.optimized.2)  # 그런데 이렇게 모형을 축소해도 다시 p-value 가 기준 이상으로 나오는 변수가 생김. 이건 설명하기가 좀 애매한데.. 공선성을 보이던 변수 하나를 제거해도 여전히 남아 있는 공선성이 이번엔 남아있는 변수에 들러붙은 것처럼 나타난다.. 이런 식? 결국 독립된 변수라는 가정이 쉽지 않은 가정이다..  

# Plot and environmental variables

both.glm.fitted <- fitted(both.glm.optimized.2) # 선형 모형과 

both.rf.fitted <- rf.both$predicted # RF 로부터 fitted 값을 추출 

origianl <- data.full$kizuki # 원자료 

both.total <- data.frame(origianl, both.rf.fitted, both.glm.fitted)

plot(both.total) # 원자료와 선형모형, 원자료와 RF 결과를 플롯 

cor(both.total) # 상관계수 행렬. 17개 그대로 쓰면 사실 여기서 GLM 결과가 더 좋아서 혼동될 수 있다. 1) GLM에서 stepwise selection 했어도 완전하게 변수 선택이 안돼서 (p-value 기준으로 허용안되는 애들이 살아남음), 그래서 over-fitted 모형이고, 그래서 R2가 무척 높다. 선형 모형의 특징인데 독립 변수 개수를 늘릴 수록 R2 는 1로 수렴함. 여튼 재밌는 지점임. 몹귀해서 안해봤는데 변수 몇 개 더 빼고 정상적인 GLM을 돌리면 적합도가 많이 떨어질 것임..  # 24일 추가. P-value 기준으로 한 번 더 선택. 여전히 GLM 결과가 더 좋게 나오긴 하나 거의 비슷해 짐. 공선성이 있는 변수를 그대로 놔두면 선형 모형에선 무슨 일이 생기는가? 추정치의 분산이 커진다. 필요한 변수가 모자라게 되면 빠진 변수의 문제가 되고, 추정치에 편의 (bias) 가 생기고 공선성 있는 변수가 남아있으면 과한 변수의 문제가 되어 추정치의 분산이 커진다. 약간 현상적으론 다를 수 있으나 대체적으로 선형 모형에선 이런 식이라 생각하면 된다. 증명은 OLS에 대한 것만 확인했음 참고로. 

# cor(both.total)

#                 origianl both.rf.fitted both.glm.fitted

# origianl        1.0000000      0.7080985       0.7173883

# both.rf.fitted  0.7080985      1.0000000       0.7720756

# both.glm.fitted 0.7173883      0.7720756       1.0000000

# Only plot data

onlyPlot.rf.fitted <- rf.onlyPlot$predicted # 사이트 측정 변수만 이용한 모형. 가장 결과가 좋았다. 대략 60% 정도의 분산을 설명했음. Pseudo-R2로 보면, cross-validated R2가 0.6 정도로 굉장히 괜찮은 결과임. 보통 RF 논문에선 이 CV R2가 0.9 넘고 하는데, 상당히 다듬어진 결과이거나 데이터 숫자가 많아서 그렇다. 

plot(origianl, onlyPlot.rf.fitted, pch=4)

## Partial dependence (or alike) # Clutter 논문에서 보여준 Partial dependence 플롯을 흉내낸 것. 제대로 하려면 데이터를 여러개 만들어내서 돌리고 합산하고 해야 하는데 여기선 단순하게 우리가 가진 predictor 중 하나에 대해 fitted 값과 산포도 그려서, 대충만 봄. 그래도 눈으로 보는 데는 괜찮음.. 논문 낼 때는 더 예쁘게.. 

# Both 

plot(x = data.both$dPC_PCA, y = both.rf.fitted, pch=5)  # 연결성에 대해 대충 선형 관계를 보임. 그런데 연결성이 2를 넘으면 오히려 감소하는 것으로 해석할 수 있는 점이 몇 개 보인다. 최소한 plateau 가 있다고 얘기해야 할 듯.  

plot(x = data.both$LogArea, y = both.rf.fitted, pch=4) # log area에 대해선 exponential 한 관계로 보이나, area가 커도 종 풍부도가 떨어지는 곳도 꽤 있다고 해석하는게 더 타당할 것 같다. 패치 크기에 로그를 취했는데도 결과와 선형 관계로는 보이지 않음, 중요 체크. 역시 RF 좋다. 

# Only plot

plot(x = data.both$siteNDVI, y = onlyPlot.rf.fitted, pch=4) # 

# 3d plot

require(scatterplot3d)

scatterplot3d(x =data.both$LogArea , y=data.both$dPC_PCA, z= both.rf.fitted, highlight.3d=T, box = FALSE, type = "h", angle = 85, pch=4)

### Proximity # Proximity 매트릭스 이용해서는 각 관측치들의 상대적인 가까움을 표시할 수 있다. 

rf.both.prox <- rf.both$proximity # 인접성은 RF  돌릴 때 부산물로 나옴. Decision tree를 만들 때 얼마나 노드들이 서로 비슷하게 분류된 횟수가 많은지를 바탕으로 함. 

rf.both.prox.mds <- cmdscale(1-rf.both.prox) # 스케일링 

# 농업/도시 지역 + 활엽수/혼합림 두 가지 정보를 동시에 출력. 명확한 효과는 보이지 않고, 이는 Variable importance 에서 두 가지 요인이 다 별로 중요치 않은 것으로 나온 것을 재확인함. 

plot(rf.both.prox.mds, col = c("green","purple")[(data.both$Matrix_type)], pch = c(1,16)[(data.both$Vegetype)], xlab="", ylab="", asp=1)

legend("bottomright", legend=c("Agri/Deci.", "Urban/Deci.", "Agi/Mixed", "Urban/Mixed"), pch=c(1,16,1,16), col=c("green", "green", "purple", "purple"))

# 도시화 정도와 서식지 타입. 범례 그려야 하는데 몹귀해서.. 서식지 타입 효과는 분명히 나타남. 

plot(rf.both.prox.mds, col = c("green","yellow", "orange", "red")[(data.both$Urbanization)], pch = c(1,16)[(data.both$Habitattype)], xlab="", ylab="", asp=1)

# legend("bottomright", legend=c("Agri/Deci.", "Urban/Deci.", "Agi/Mixed", "Urban/Mixed"), pch=c(1,16,1,16), col=c("green", "green", "purple", "purple"))

# 반응 표면 그리기 위해 데이터 준비. 

data3d <- data.frame(x =data.both$LogArea , y=data.both$dPC_PCA, z= both.rf.fitted)

# data3d