R 프로그램을 이용해서 다중공선성 확인하기 - 2편

이번 글은 이전 글(2023.03.15 - [머신러닝 & 딥러닝] - R기반 래스터 다중공선성 확인하기)에 이어서 래스터 자료를 이용해서 다중공선성을 확인하는 방법 중 생물종 위치나 산불, 산사태, 범죄 등 특정 사건이 발생한 위치에 해당하는 독립변수 값만을 이용해서 한다. 또한 상관계수가 아니 분산팽창계수(VIF)를 이용한다.

 

일단 작업폴더를 설정하고 필요한 라이브러리를 설치하고 불러온다

# 작업폴더 설정
setwd("c:/R_work")
getwd()

# 라이브러리 불러오기
library(dplyr)
library(raster)

남색이마잠자리 위치에서 생물기후변수 19개 사이의 상호관계를 확인하여 다중공선성을 검토한다.

이를 위해서 남색이마잠자리 위치관련 shp파일을 불러오고 생물기후변수 19개 변수를 불러온다.

#남색이마잠자리 shp 파일 지정
file <- "./남색이마잠자리_pt.shp

# shape 파일 불러오기
sp <- shapefile(file) 
head(sp)
class(sp)

# bioclim data 해상도 10분짜리 다운로드하기

bio <- raster::getData('worldclim', var='bio', res=10)
bio
names(bio)

# bio 변수 1번인 평균기온 출력하고 그 위에 남색이마잠자리 위치포인터 올리기
plot(bio[[1]])  
points(sp)

동남아시아로 한정하기 위해서 대상지역에 해당하는 shp 파일을 불러오고 대상지역에 x좌표와 y좌표의 최소값과 최대값을 이용한다.

ext <- "./남색이마잠자리_AOI.shp" # shape 파일 지정

sp_ext <- shapefile(ext) # shape 파일 불러오기
head(sp_ext)
class(sp_ext)

e <- extent(sp_ext)
e

대상지역이 결정되면 대상지역 안에 있는 생물종위치와 환경자료를 잘라낸다.

spg <- crop(sp, e)
points(spg,col='red')
bioc <- crop(bio, e)
plot(bioc[[1]])
points(spg)
class(spg)

래스터 파일로 구성된 환경변수 자료의 상호 상관관계를 계산하기 위해서 raster.cor.matrix()함수를 이용한다.

또한 분산팽창계수인 VIF지수를 사용하기 위해서는 usdm 라이브러리를 불러오고 vif() 함수를 이용한다.

다음은 래스터 전체를 대상으로 상관관계와 분산팽계수를 계산한다.

# 다중공선성 영향 제거
library(usdm)
raster.cor.matrix(bioc)

#------------------------
vif(bioc)
#------------------------------

다음 VIF 지수값을 보면 10이상의 값을 가진 환경변수를 볼 수 있다. VIF 지수값이 10이상이면 다중공선성이 있다고 할 수 있다.

이제 생물종의 위치에서 래스터 환경변수 값을 raster::extract()함수를 이용하여 추출한 다음 상관관계와 분산팽창계수를 계산한다.

ex <- raster::extract(bioc,spg)
class(ex)
head(ex)
ex
nrow(ex)

#-------------------------
v <- vifstep(ex)
v

VIF 계산한 결과 지수값이 10이상 독립변수를 제외한 다음 8개의 독립변수가 선택된 것을 알 수 있다

 

위에서 분산팽창계수 VIF 값이 10이하인 독립변수인 래스터 파일만을 exclude()함수를 이용해서 선택한다.

bioc_ex_vif <- exclude(bioc, v)
plot(bioc_ex_vif)

다음은 참고로 독립변수 간의 상관관계를 계산하는 방법이다. matrxi array 자료 형태인 ex 파일을 data.frame자료 형태로 변환한 다음 cor()함수를 이용해서 상관계수를 계산한다.

#-----------------------
ex_df =as.data.frame(bioc_ex_vif)
ex_df
class(ex_df)
cor(ex_df, use = "pairwise")

다음은 QGIS이나 다른 목적으로 사용하기 위해서 다중공선성이 없고 생물종 위치에 해당하는 독립변수 값을 추출하여 GIS shp파일을 만드는 구문이다.

#--------------
ex_sp <- raster::extract(bioc_ex_vif,spg, sp=TRUE)
head(ex_sp)
raster::shapefile(ex_sp, "ex_val.shp", overwrite=TRUE)

전체코드

# 작업폴더 설정
setwd("c:/R_work")
getwd()


# 라이브러리 불러오기

library(dplyr)
library(raster)


#남색이마잠자리 shp 파일 지정
file <- "./남색이마잠자리_pt.shp" 

# shape 파일 불러오기
sp <- shapefile(file) 
head(sp)
class(sp)


# bioclim data 해상도 10분짜리 다운로드하기

bio <- raster::getData('worldclim', var='bio', res=10)
bio
names(bio)

# bio 변수 1번인 평균기온 출력하고 그 위에 남색이마잠자리 위치포인터 올리기
plot(bio[[1]])  
points(sp)


class(bio[[1]])

ext <- "./남색이마잠자리_AOI.shp" # shape 파일 지정

sp_ext <- shapefile(ext) # shape 파일 불러오기
head(sp_ext)
class(sp_ext)

e <- extent(sp_ext)
e

spg <- crop(sp, e)
points(spg,col='red')
bioc <- crop(bio, e)
plot(bioc[[1]])
points(spg)
class(spg)

for (x in 1:19) {
  print(names(bioc[[x]]))
  writeRaster(bioc[[x]], paste('wc10_asc/','bio_',x,'.asc'), NAflag=-9999, overwrite=TRUE)
}

# 다중공선성 영향 제거
library(usdm)
raster.cor.matrix(bioc)

#------------------------
vif(bioc)
#------------------------------

ex <- raster::extract(bioc,spg)
class(ex)
head(ex)
ex
nrow(ex)

#-------------------------
v <- vifstep(ex)
#vifcor
v
bioc_ex_vif <- exclude(bioc, v)
plot(bioc_ex_vif)

#-----------------------
ex_df =as.data.frame(ex)
ex_df
class(ex_df)
cor(ex_df, use = "pairwise")

#--------------
ex_sp <- raster::extract(bioc_ex_vif,spg, sp=TRUE)
head(ex_sp)
raster::shapefile(ex_sp, "ex_val.shp", overwrite=TRUE)