# HSM(Habitat Suitability Model) 2021/02/11
## 次回
https://hackmd.io/lenjCgv1S7-HKiJEB_GIeA
## Tasks by 2/11
- 6章：中嶋(-6.2.6)、向峯(6.2.7-6.3.2)、八尾(6.3.3-)
- 7章：岡本(7の4-)、八木(7の1-3)
を読んできて資料を作る。ソースコードを走らせる。
各章の中での配分は各章ごとに決めましょう。
6章は
6.0.0 - 6.2.6：中嶋
6.2.7 - 6.3.2：向峯
6.3.3 - 6.4.3：八尾
とかでどうでしょう？
## 6.0
"raster"パッケージは探索的データ解析（explanatory analysis）や空間モデリングに適している．
> 探索的データ分析(EDA)は，要約統計量およびグラフ化ツールを使用して、データを把握し，データから得られる知見を理解する調査プロセスである． (https://www.jmp.com/ja_jp/online-statistics-course/exploratory-data-analysis.html)
## 6.1
#### 数値標高モデル (陸域表面を3次元表示する)
1.[ GLOBE (GTOPO30)](https://www.ngdc.noaa.gov/mgg/topo/gltiles.html)
地理的投影で30秒（赤道付近で1 km）の解像度
全球をカバーするが高精度ではない
2. [ETOPO1](https://www.ngdc.noaa.gov/mgg/global/global.html)
地理的投影で1分（赤道付近で2 km）の解像度
水深を含むが，GLOBEと同様のデータ
3. [SRTM](https://srtm.csi.cgiar.org/)
全球で利用できる90 mの解像度のデータセット
60°S-60°Nまで途切れなしで利用できる
4.[ASTER-DEM](https://asterweb.jpl.nasa.gov/gdem.asp)
空間解像度は30 mで1°のグリッドセルで全球レベルで入手できる
83°S-83°Nまで途切れなしで利用できる
#### 気象データ
1. [WorldClim](https://worldclim.org/)
地理的投影で30秒(約1 km)の空間解像度で，粗い精度でも（2.5分，5分，10分）が利用可能
* 月平均降水量
* 最高気温
* 最低気温
---
(粗い精度2.5分，5分，10分は上記＋)
* 平均気温
* solar radiation
* wind speed
* water vapour pressure
上に加えて[生物気候変数(bioclim)](https://worldclim.org/data/bioclim.html)も提供する．
2. [PRISM](https://prism.oregonstate.edu/)
1985-2010年で，4 kmの空間解像度で月と年で利用できる
1981-2010年の気候学的平年値は800 mと4 kmの解像度で利用できる
1981年からの日データも4kmの解像度で提供している
3. [DAYMET](https://daymet.ornl.gov/)
1980-1997年の最高最低気温，降水量，全天日射短波放射，飽差について毎日の地表の気象をマッピングする
NetCDFフォーマットで1kmの空間解像度で1980-2015年の期間をカバーする
北米，プエルトリコ，ハワイを網羅する
上記全ては1 kmの解像度で入手可能（より粗いデータだとより多くのデータがある）
### 土地被覆データ
土地被覆及び植生被覆のデータセット
1. 欧州宇宙機関(ESA)
地球規模の土地被覆データ
2. National land cover database
アメリカの土地被覆
主題分類，不透水度，樹冠被覆率など
3． Coordinated information on the europe environment
ヨーロッパ圏
土地被覆．土地利用
### 国境政治単位
1. GADM
国などの世界の行政区域及び州・県などの空間データ
2. Natural Earth
国，紛争地域，県，州などの一次区分
人口密集地，都市ポリゴン，公園および保護地域，水域境界
3. ESPON
ヨーロッパ
4. Marine region
## 6.2
重いのか，plotするとRが頻繁に反応しなくなる
```
par(mfrow=c(3,1))
plot(elev1,col=rev(topo.colors(50)),main="Elevation")
plot(bio3,col=heat.colors(100),main="Bio.3")
plot(bio11,col=rainbow(100),main="Bio.11")
par(mfrow=c(1,1))
```
```
par(mfrow=c(2,2))
plot(bio3,bio12,xlab="bio3",ylab="bio12",col="gray55")
plot(bio3,bio7,xlab="bio3",ylab="bio7",col="gray55")
plot(bio3,bio11,xlab="bio3",ylab="bio11",col="gray55")
plot(bio3,elev1,xlab="bio3",ylab="elevation",col="gray55")
par(mfrow=c(1,1))
```
```
plot(elev1, main="Elevation Contours")
contour(elev, nlevels=7, levels=c(0, 500, 1000, 1500, 2000, 3000, 4000, 5000), add=TRUE, labels="", lwd=.2)
```
```
elev_sa <- crop(elev, extent(-85,-30,-60,15))
elev_na <- crop(elev, extent(-188,-50,15,90))
plot(elev_sa, main="Elevation Contours South America ")
contour(elev_sa, nlevels=7, levels=c(0, 500, 1000, 1500, 2000, 3000, 4000, 5000), add=TRUE, labels="", lwd=.2)
```
```
plot(elev_na,col=terrain.colors(40),main="Elevation Contours")
lines(iso1, lwd=0.2)
```
```
lat<-raster("raster/other/latitude.tif")
lon<-raster("raster/other/longitude.tif")
tcold<-4.138e+02 + (-3.624e-02 * elev1) + (-8.216e+00 * abs(lat)) +
(-1.794e+00 * abs(lon)) + (7.122e-03 * lon^2)
diff_obs_model_temp <-bio11 - tcold
par(mfrow=c(2,1))
plot(tcold, col=rev(rainbow(100))[20:100],main="Modelled mean
temperature of the coldest quarter")
contour(elev, nlevels=7, levels=c(0, 500, 1000, 1500, 2000, 3000, 4000,
5000), add=T, labels="", lwd=.3)
plot(diff_obs_model_temp, col=rev(rainbow(100))[20:100],main=
"Difference between modelled and observed temperatures")
contour(elev, nlevels=7, levels=c(0, 500, 1000, 1500, 2000, 3000, 4000,
5000), add=T, labels="", lwd=.3)
par(mfrow=c(1,1))
```
### 6.2.7
グリッドは空間的に重ねる（スタックする）ことができる。
スタックしたグリッドに対しての処理は、全てのグリッドへ一括で処理を行うことができる。
`stack()`
`aggreagate()`
### 6.2.8
測地系が異なると投影時にズレが生じる。->誤った解析をしてしまう可能性も。
なので処理を行う前に測地系を揃える必要がある。
関数は`raster::projectRaster()`
### 6.2.9
種分布データがデータフレームである場合、スタックするためにspatialクラスに設定する必要がある。該当部分は以下。
```r=
pinus_edulis <- read.table("tabular/species/pinus_edulis_occ.csv", sep=",",header=TRUE)
coordinates(pinus_edulis) <- c("lon", "lat")
projection(pinus_edulis) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
```
1行目でデータフレームの読み込み
2行目でspatialクラスへの変換
3行目で座標系の指定
を行っている。
種のデータはgbifからダウンロードして使うことも可能。
{dismo}を使う（{rgbif}でも可能っぽい）。
```r=
pinus_edulis <- gbif('pinus', 'edulis', download=T, geo=T, sp=T, removeZeros=T)
names(pinus_edulis)[1]<-"dwnld.date"
```
### 6.2.10
GLMまでをお試しでやってみる。
イメージとしては重ね合わせたグリッドを説明変数/応答変数に割り当てて、地点数分のサンプルがある感じ。
```r=
dat <- tribble(~lat_lon, ~distribution, ~temperature, ~precipitation,
1, 0, 24.5, 40,
2, 0, 25.1, 60,
3, 1, 30.0, 50,
4, 1, 31.2, 40)
glm(distribution ~ temperature + precipitation, data=dat) %>% summary
> Call:
> glm(formula = distribution ~ temperature + precipitation, data = dat)
> Deviance Residuals:
> 1 2 3 4
> 0.03205 -0.05804 0.11608 -0.09009
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) -4.1146405 1.0502104 -3.918 0.159
> temperature 0.1674848 0.0287320 5.829 0.108
> precipitation -0.0005198 0.0101826 -0.051 0.968
> (Dispersion parameter for gaussian family taken to be 0.02598903)
> Null deviance: 1.000000 on 3 degrees of freedom
> Residual deviance: 0.025989 on 1 degrees of freedom
> AIC: -0.79399
>
> Number of Fisher Scoring iterations: 2
```
- ステップワイズ変数選択：モデルを少しずつ変えながらAICなりの基準で最適なモデルを探す方法。フルモデルから始める場合、forward-backward（変数増減法）、ヌルモデルから始める場合backward-forward（変数減増法）。
- adjusted D-2：自由度によって調整された分散の値。小さい程よいモデル。AICと同じでこの値だけで判断するのはよくない。
### 6.3.1
リモセンは有用。
### 6.3.2
Landsatの説明。
Landsat7とLandsat8とでバンドが異なる。
<details>
<summary>tips (OS X)</summary>
6章を走らせた時のエラー集。
```r=
library(maptools)
Checking rgeos availability: FALSE
Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
which has a restricted licence. It is disabled by default;
to enable gpclib, type gpclibPermit()
```
なんかエラー出たので
https://stackoverflow.com/questions/21093399/how-to-turn-gpclibpermit-to-true
に従って
```r=
install.packages("gpclib")
library(gpclib)
library(maptools)
```
これでok.
```r=
elev <- raster("raster/topo/GTOPO30.tif")
Error in .local(.Object, ...) :
.rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", でエラー:
Cannot create a RasterLayer object from this file.
```
エラーが出たので本の脚注に従って
https://github.com/hsdm-hub/hsdm/blob/master/data/raster/topo/GTOPO30
からデータをダウンロードし直した。これでok.
```r=
pts.clim<-raster::extract(world.stk, pinus_edulis, method="bilinear")
.xyValues(x, as.matrix(y), ...) でエラー: xy should have 2 columns only.
```
extract()がconflictしているらしい。めんどくさいがraster::extract()で対応。conflicted::conflict_preferを使っても良い。
```r=
install.packages("ecospat")
```
でエラー。何やら{deldir}をインストールするためにgfortranがないらしい。下記の記事でも生じている。
https://stackoverflow.com/questions/25315422/cant-find-gfortran-4-8-to-build-package
```bash=
brew install gfortran
```
で対応。
ただそれでも
```r=
install.packages("deldir")
...
...
ld: warning: directory not found for option '-L/usr/local/gfortran/lib/gcc/x86_64-apple-darwin18/8.2.0'
ld: warning: directory not found for option '-L/usr/local/gfortran/lib'
ld: library not found for -lgfortran
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make: *** [deldir.so] Error 1
ERROR: compilation failed for package ‘deldir’
* removing ‘/Library/Frameworks/R.framework/Versions/4.0/Resources/library/deldir’
```
となる。
解決方法は{deldir}のバイナリ版を選べば良い（盲点だった...）。これは
```r=
install.packages("deldir")
```
の際に"Yes"ではなく"no"を選択すれば良い。
</details>
## 6.3
### 6.3.3
Landsatのデータ取得
米国地質調査所USGSのHPからLandsatの衛星写真を取得できる．無料の利用登録が必要．
↓参考HP
- [オンラインデータベースからの衛星画像の取得](https://learn.arcgis.com/ja/projects/download-imagery-from-an-online-database/)
- [Landsat 8 画像を使ってみよう](https://blog.esrij.com/2017/04/28/post-24474/)
p.95に注意書きとして，「Landsatのデータが更新されており，以下のデータは取得できない．」とあったが，USGSのHPから取得できた．
→取得したデータをlandsat_self2というフォルダに保存した．
- スイス チューリッヒ州を包括するシェープファイルを読み込んで画像をトリミング
```
Cantons <- readShapePoly("vector/swiss/Swiss_Cantons.shp")
警告メッセージ:
readShapePoly is deprecated; use rgdal::readOGR or sf::st_read
Cantons <- rgdal::readOGR("vector/swiss/Swiss_Cantons.shp")
OGR data source with driver: ESRI Shapefile
Source: "/Users/yaoakifumi/Desktop/informatics/hsm/Ch.6/vector/swiss/Swiss_Cantons.shp", layer: "Swiss_Cantons"
with 48 features
It has 8 fields
Zurich <- Cantons[Cantons$NAME == "ZUERICH",]
```
maptool::readShapePolyだとエラーが出たので，rgdal::readOGRを使用した．
- Landsatバンドの読み込み
```
#load data
band1_blue <- raster("landsat_self2/LE07_L1TP_194027_20010824_20170203_01_T1_B1.TIF")
band2_green <- raster("landsat_self2/LE07_L1TP_194027_20010824_20170203_01_T1_B2.TIF")
band3_red <- raster("landsat_self2/LE07_L1TP_194027_20010824_20170203_01_T1_B3.TIF")
band4_nir <- raster("landsat_self2/LE07_L1TP_194027_20010824_20170203_01_T1_B4.TIF")
band5_swir1 <- raster("landsat_self2/LE07_L1TP_194027_20010824_20170203_01_T1_B5.TIF")
band7_swir2 <- raster("landsat_self2/LE07_L1TP_194027_20010824_20170203_01_T1_B7.TIF")
```
- Landsatバンドをチューリッヒ州に合わせてトリミング
```
#triming
band1_blue_crop <- crop(band1_blue, extent(Zurich))
band2_green_crop <- crop(band2_green, extent(Zurich))
band3_red_crop <- crop(band3_red, extent(Zurich))
band4_nir_crop <- crop(band4_nir, extent(Zurich))
band5_swir1_crop <- crop(band5_swir1, extent(Zurich))
band7_swir2_crop <- crop(band7_swir2, extent(Zurich))
```
- 6つの画像を一つの複合レイヤーに結合
```
#connect layers
L7 <- brick(band1_blue_crop, band2_green_crop, band3_red_crop,
band4_nir_crop, band5_swir1_crop, band7_swir2_crop)
```
- トリミング〜複合レイヤーへの積み上げは以下のコードでもできる
```
#more effective code
tmp <- stack(band1_blue, band2_green, band3_red,
band4_nir, band5_swir1, band7_swir2)
L7 <- crop(tmp, extent(Zurich))
```
- レイヤー名の割り当て
```
#apply layer names
names(L7) <- c("band1_blue", "band2_green", "band3_red",
"band4_nir", "band5_swir1", "band7_swir2")
names(L7)
```
- カラー合成画像
```
#plot color
par(mfrow = c(1,3))
plotRGB(L7, 3, 2, 1, stretch = "lin")
plotRGB(L7, 4, 3, 2, stretch = "lin")
plotRGB(L7, 6, 4, 3, stretch = "lin")
par(mfrow = c(1,1))
```
### 6.3.4 生態学的分析のためのデータ処理
リモートセンシングでよく用いられる指数をバンド画像から再計算する．
- 指数
リモートセンシングで指数を利用するのは，衛星画像のバンドそのままでは見えない特徴を強化して分析可能にするためである．
代表的な指数としては
・<b>正規化植生指数NDVI
・単純比SR
・土壌調整植生指数SAVI
・正規化水指標NDWI</b>
がある．
また，植生と都市開発の変化を分析するためのタッセルドキャップ変換指標がある．
- 代表的な4指標の計算
```
NDVI <- (L7$band4_nir - L7$band3_red) / (
L7$band4_nir + L7$band3_red)
NDWI <- (L7$band4_nir - L7$band5_swir1) / (
L7$band4_nir + L7$band5_swir1)
SR <- L7$band4_nir / L7$band3_red
```
SAVIと後述するタッセルドキャップの計算には，At-Sensor反射率または表面反射率が必要なため，landsatパッケージを用いて算出する．
landsatパッケージに対応できるようにSpatialGridDataFrameにデータ変換
```
#format change
L7_sp <- as(L7, "SpatialGridDataFrame")
L7_sp1 <- L7_sp[1]
L7_sp2 <- L7_sp[2]
L7_sp3 <- L7_sp[3]
L7_sp4 <- L7_sp[4]
L7_sp5 <- L7_sp[5]
L7_sp7 <- L7_sp[6]
```
landsat::radiocorrでAt-Sensor反射率を計算
```
#package
library(landsat)
#calculate At-Sensor
L7_ref_sp1 <- radiocorr(L7_sp1, Grescale = 0.76282,
Brescale =-1.52, sunelev = 48.29,
edist = ESdist("2011-08-24"),
Esun = 1957, method = "apparentreflectance")
L7_ref_sp2 <- radiocorr(L7_sp2, Grescale = 1.44251,
Brescale =-2.84, sunelev = 48.29,
edist = ESdist("2011-08-24"),
Esun = 1826, method = "apparentreflectance")
L7_ref_sp3 <- radiocorr(L7_sp3, Grescale = 1.03988,
Brescale =-1.17, sunelev = 48.29,
edist = ESdist("2011-08-24"),
Esun = 1554, method = "apparentreflectance")
L7_ref_sp4 <- radiocorr(L7_sp4, Grescale = 0.87258,
Brescale =-1.51, sunelev = 48.29,
edist = ESdist("2011-08-24"),
Esun = 1036, method = "apparentreflectance")
L7_ref_sp5 <- radiocorr(L7_sp5, Grescale = 0.11988,
Brescale =-0.37, sunelev = 48.29,
edist = ESdist("2011-08-24"),
Esun = 215, method = "apparentreflectance")
L7_ref_sp7 <- radiocorr(L7_sp7, Grescale = 0.06529,
Brescale =-0.15, sunelev = 48.29,
edist = ESdist("2011-08-24"),
Esun = 80.67, method = "apparentreflectance")
```
SAVI
緑の植生被覆の量に応じて変化するパラメータLが必要であり，高密度植生地域ではL = 0, 緑の植生のない地域ではL = 1を使用する．
デフォルトのL = 0.5で大抵うまくいく．
```
#SAVI
L <- 0.5
SAVI <- ((raster(L7_refl_sp4) -
raster(L7_refl_sp3))/
(raster(L7_refl_sp4) +
raster(L7_refl_sp3) + L) ) * (1 + L)
```
- タッセルドキャップ変換
landsat::tasscapで計算する．
```
#tasscap
L7_tc <- tasscap("L7_refl_sp", sat = 7)
L7_Brightness <- raster(L7_tc[[1]])
L7_Greenness <- raster(L7_tc[[2]])
L7_Wetness <- raster(L7_tc[[3]])
```
生のバンドよりも生態学的に意味のある<b>明度，緑，湿り度</b>の3つのレイヤーを抽出できる．
- 分析結果の描画
250 mの解像度の陰影を利用して，分析結果を描画する．
主要な4つの指数
```
par(mfcol = c(2,2))
plot(hill_250m_utm, col = grey(0:100/100),
legend = F, axes = F,
ext = extent(SAVI), main = "NAVI")
plot(NDVI, col = rev(terrain.colors(20, alpha = 0.6)), add = T)
plot(hill_250m_utm, col = grey(0:100/100),
legend = F, axes = F,
ext = extent(SAVI), main = "SR")
plot(SR, col = rev(terrain.colors(20, alpha = 0.6)), add = T)
plot(hill_250m_utm, col = grey(0:100/100),
legend = F, axes = F,
ext = extent(SAVI), main = "NDWI")
plot(NDWI, col = rev(topo.colors(20, alpha = 0.6)), add = T)
plot(hill_250m_utm, col = grey(0:100/100),
legend = F, axes = F,
ext = extent(SAVI), main = "SAVI")
plot(SAVI, col = rev(terrain.colors(20, alpha = 0.6)), add = T)
```
タッセルドキャップ変換
```
par(mfrow = c(1,3))
plot(hill_250m_utm, col = grey(0:100/100), legend = F, axes = F,
ext = extent(SAVI), main = "Brightness")
plot(L7_Brightness, col = rev(paste(ygb.c(20), "B3",
sep = "")), add = T)
plot(hill_250m_utm, col = grey(0:100/100), legned = F, axes = F,
ext = extent(SAVI), main = "Greenness")
plot(L7_Greenness, col = rev(terrain.colors(20, alpha = 0.6)),
add = T)
plot(hill_250m_utm, col = grey(0:100/100), legend = F, axes = F,
ext = extent(SAVI), main = "Wetness")
plot(L7_Wetness, col = rev(topo.colors(20, alpha = 0.6))
, add = T)
```
## 6.4
### 6.4.1
・予測変数は対象種の分布と因果関係がなければならない．
・相対湿度や霜日の頻度などといった予測変数は，地形地図から導出されることが多いが，誤差伝播によりシンプルな地形地図よりも誤差が蓄積することに留意する必要がある．
・間接的な要因を予測変数とすることは，原則として避けるべきである．(e.g. 森林限界の予測に気温ではなく標高を用いるなど)
### 6.4.2
・予測変数間に強い相関がある場合，多重共線性が問題となる場合が多い．
・この問題を回避するために，変数間の相関をまず可視化する．
```
library(ecospat)
data <- read.csv("tabular/bioclim/current/bioclim_table.csv",
header = TRUE, sep = ",")
ecospat.cor.plot(data[,4:8])
```
可視化すると，
・bio4 (気候の季節性) とbio7 (年間気候幅) が強く相関する (r = 0.97) .
・bio4 (気候の季節性) とbio11 (最も寒い四半期の平均気温) も強く相関する (r = 0.93) ．
・bio7は正規分布に近いが，bio12は歪んだ分布を示す．
ことなどがわかる．
多重共線性の問題を緩和するには，強く相関する変数を取り除くことが有効である．
→ 一般にr = 0.8を基準とする．
ただし，ペアワイズ相関図では見えない相関構造がある場合があるため，分布拡大係数 (variance inflation factor ,VIF) を使う．
VIFは単なるペアワイズな相関関係ではなく，変数間の線形相関構造を検出できる．
VIFの使用例
```
library(usdm)
vif(data[,4:8])
Variables VIF
1 bio3 6.780899
2 bio4 61.448174
3 bio7 31.811787
4 bio11 11.606773
5 bio12 2.111669
```
5~10を超えるVIF値は致命的な多重共線性があることを示す．ここでは，bio4が高いVIF値を示す．
そこで，bio4を取り除いてみる．
```
vif(data[,c(4, 6:8)])
Variables VIF
1 bio3 5.809681
2 bio7 4.662858
3 bio11 6.902710
4 bio12 1.948251
```
全ての変数のVIF値が6程度かそれ以下となり，解析に使用できる．
usdmパッケージのvifstep関数を使用すると同じような結果を得ることができる．
```
# command in hsm
>vifstep(world.stk[[2:5]])
h(simpleError(msg, call)) でエラー:
引数 'x' の評価中にエラーが起きました (関数 'vifstep' に対するメソッドの選択時): オブジェクト 'world.stk' がありません
#my code
>vifstep(data[4:8], th = 10)
1 variables from the 5 input variables have collinearity problem:
bio4
After excluding the collinear variables, the linear correlation coefficients ranges between:
min correlation ( bio12 ~ bio11 ): 0.5014162
max correlation ( bio11 ~ bio3 ): 0.8911935
---------- VIFs of the remained variables --------
Variables VIF
1 bio3 5.751689
2 bio7 4.668901
3 bio11 6.820122
4 bio12 1.908185
```
あるレベルの相関を許容したとき (e.g. r = 0.7) に，どの変数が残るかをvifcorで確認できる．
```
vifcor(data[,4:8], th = .7)
3 variables from the 5 input variables have collinearity problem:
bio4 bio11 bio7
After excluding the collinear variables, the linear correlation coefficients ranges between:
min correlation ( bio12 ~ bio3 ): 0.6033323
max correlation ( bio12 ~ bio3 ): 0.6033323
---------- VIFs of the remained variables --------
Variables VIF
1 bio3 1.572351
2 bio12 1.572351
```
予測変数間の相関に注意する理由は
・変数間に相関があると，変数の影響を正確に推定できないため
・変数間の相関構造は変化する可能性があるが，強い相関のある変数をモデルに入れすぎると，シミュレーション結果がその相関によって制限され，誤った予測をしてしまう恐れがあるため
### 6.4.3
豊富にある空間データから，モデリングに使用する変数を選択する必要がある．
・ただし，変数の事前選択法については未解決な麺が多い．
・一般的には，重要でない変数の重みを下げたり除いたり，強く相関する変数を除くことが行われる．
よく行われる手順2つ
・最も重要な変数 (群) を選択し，互いに強すぎる相関がなくVIFによる検定をパスする説明変数を選択する．
・概念的な観点や実際の実験結果等の生態学的な知見から分布に影響しうる説明変数を選択する．この場合も変数間の相関には注意する必要がある．
あまり注目されていないが，変数変換も重要な側面である．
```
cor(data$bio12, data$bio7)
[1] -0.6439925
```
```
cor(log(data$bio12 + 0.001), data$bio7)
[1] -0.4410228
```
左に歪んだ分布を示すbio12を対数変換すると，bio12とbio7の間の相関関係は小さくなる．
## 7.1-7.3
担当:八木
https://rpubs.com/h8gi/727390
## 7.4
担当:岡本
https://rpubs.com/okam/724338
Sign in with Wallet
Connect another wallet