R で possum-foraging の基礎分析

# はじめに
[ポッサム](https://ja.wikipedia.org/wiki/%E3%83%9D%E3%83%83%E3%82%B5%E3%83%A0)というニュージーランドなどに生息する動物にGPSセンサーを付け、行動記録を取ったデータセット（[possum-foraging](https://data.world/nz-govt/possum-foraging)）がある。この基礎分析を行ったので、結果とともに基礎分析で分かった事、そこからデータに行った修正について記載していく。

# パッケージとデータの読み込み

```R:コード
library(dplyr)
library(ggplot2)

df <- read.csv("./yockney-et-al-data-for-poned1439135.csv")
str(df)
```

```:出力
'data.frame':	33581 obs. of  8 variables:
 $ TrackID      : int  1 2 3 4 5 6 7 8 9 10 ...
 $ Study.Site   : chr  "NC" "NC" "NC" "NC" ...
 $ Animal.code  : chr  "A8301" "A8301" "A8301" "A8301" ...
 $ NZDT_datetime: chr  "21/09/2009 0:56" "21/09/2009 2:57" "21/09/2009 4:57" "21/09/2009 7:56" ...
 $ Latitude     : num  -42.1 -42.1 -42.1 -42.1 -42.1 ...
 $ Longitude    : num  174 174 174 174 174 ...
 $ HDOP         : num  1.8 9.1 2.1 1.5 4.9 2.9 3 2.5 4.1 8 ...
 $ Season       : chr  "spring" "spring" "spring" "spring" ...
```

いくつかの文字列と位置情報がふくまれたデータであることが確認できた。


# 量的変数の統計量の算出
これまでの分析と同様にまずデータ項目の基礎統計量を算出する。

```R:
classes <- df %>% summarise_all(class) %>% t()
means <- df %>% select(!where(is.character)) %>% summarise_all(mean, na.rm=TRUE) %>% t()
vars <- df %>% select(!where(is.character)) %>% summarise_all(var, na.rm=TRUE) %>% t()
sds <- df %>% select(!where(is.character)) %>% summarise_all(sd, na.rm=TRUE) %>% t()
qs <- df %>% select(!where(is.character)) %>% summarise_all(quantile, na.rm=TRUE) %>% t()
colnames(qs) <- c("min", "q14", "median", "q34", "max")
is.nas <- df %>% is.na %>% as.data.frame %>% summarise_all(sum) %>% t()

stat.1 <- data.frame(
  rownames=rownames(classes),
  classes=classes,
  is.nas=is.nas,
  is.not.nas=nrow(df)-is.nas
)

stat.2 <- data.frame(
  rownames=rownames(means),
  means=means,
  vars=vars,
  sds=sds,
  qs
)

stat <- full_join(stat.1, stat.2, by="rownames")
write.csv(x=stat, file="summary.csv")
```

出力の summary.csv は以下のような出力となる

| |rownames|classes|is.nas|is.not.nas|means|vars|sds|min|q14|median|q34|max|
|:--|:--|:--|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|
|1|TrackID|integer|0|33581|16791|93976428.5|9694.144031|1|8396|16791|25186|33581|
|2	|Study.Site	|character	|0	|33581	|NA	|NA	|NA	|NA	|NA	|NA	|NA	|NA|
|3	|Animal.code	|character	|0	|33581	|NA	|NA	|NA	|NA	|NA	|NA	|NA	|NA|
|4	|NZDT_datetime	|character	|0	|33581	|NA	|NA	|NA	|NA	|NA	|NA	|NA	|NA|
|5	|Latitude	|numeric	|0	|33581	|-16.17020564	|1725.902724	|41.54398541	|-45.30308	|-45.28762	|-42.14543	|45.28499	|45.30198|
|6	|Longitude	|numeric	|0	|33581	|170.2712649	|3.213393283	|1.792594009	|169.26576	|169.28312	|169.29043	|169.30543	|173.56041|
|7	|HDOP	|numeric	|0	|33581	|3.071507698	|3.44993206	|1.857399273	|0.8	|1.9	|2.6	|3.6	|12.7|
|8	|Season	|character	|0	|33581	|NA	|NA	|NA	|NA	|NA	|NA	|NA	|NA|

ここで注目すべきは Latitude の q34, max の項目。値がプラスになっていて明らかにニュージーランドから出てしまっているのが分かる。GPSの値違いかデータの一部が±反転してしまっているよう。その値に引っ張られる形で Latitude の means もニュージーランドではなくなってしまっている。Rの`leaflet`を用いてこの means の座標をプロットしてみると、

```R:コード
library(leaflet)
m <- leaflet() %>% addTiles()
m <- m %>% 
  addCircles(lng = means[3], lat = means[2]) %>% 
  setView(lng = means[3], lat = means[2], zoom = 4)
m
```
![R_map_miss.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/cb0a5c7f-5747-4850-1700-5c63b074c7e8.png)

となり、ニュージーランドの北の海を示す事になる。陸上生物かつ羽もないポッサムの分布としてはおかしいことが分かる。

__以上より、Latitude のプラスのデータをマイナスに直すという処理が必要になる。__


# 質的変数の度数分布の確認
時刻の列を除くと、Study.Site, Animal.code, Season の三つの質的変数がある。それぞれについて度数分布を確認していく。

```R:コード
g.Study.Site <- ggplot(df, aes(x = Study.Site))
g.Study.Site <- g.Study.Site + geom_bar()
plot(g.Study.Site)
```
![R_Study_Site.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/fea797af-ca1c-b13a-1197-e83522336d1e.png)

```R:コード
g.Animal.code <- ggplot(df, aes(x = Animal.code))
g.Animal.code <- g.Animal.code + geom_bar()
g.Animal.code <- g.Animal.code + theme(axis.text.x = element_text(angle = 90, hjust = 1))
plot(g.Animal.code)
```
![R_Animal_code_miss.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/e73c8031-9dc2-f796-82be-ace3699100b1.png)

```R:コード
g.Season <- ggplot(df, aes(x = Season))
g.Season <- g.Season + geom_bar()
plot(g.Season)
```
![R_Season_miss.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/1f621609-6245-73e1-d2e3-3c49827dd869.png)

Animal.code と、Season でそれぞれおかしいところがある。Animal.code では、 kristin2 と Kristin2、Season では、spring と Spring、summer と Summer のそれぞれ小文字と大文字との表記ゆれがあった。

__これらのデータについても統一して処理する必要がある。__

# ±の変換、表記ゆれの修正
上記で見つけたデータの修正点を修正するコードを実装していく。

```R:コード
for (i in 1:nrow(df)) {
  if (df$Latitude[i] > 0) {
    df$Latitude[i] <- df$Latitude[i] * -1 #±の変換
  }
  if (df$Animal.code[i] == "kristin2") {
    df$Animal.code[i] <- "Kristin2" #表記ゆれの修正
  }
  if (df$Season[i] == "spring") {
    df$Season[i] <- "Spring" #表記ゆれの修正
  } else if (df$Season[i] == "summer") {
    df$Season[i] <- "Summer" #表記ゆれの修正
  }
}
```

修正ができているか再度、統計量を算出、プロットを作って確認する。

```R:
classes2 <- df %>% summarise_all(class) %>% t()
means2 <- df %>% select(!where(is.character)) %>% summarise_all(mean, na.rm=TRUE) %>% t()
vars2 <- df %>% select(!where(is.character)) %>% summarise_all(var, na.rm=TRUE) %>% t()
sds2 <- df %>% select(!where(is.character)) %>% summarise_all(sd, na.rm=TRUE) %>% t()
qs2 <- df %>% select(!where(is.character)) %>% summarise_all(quantile, na.rm=TRUE) %>% t()
colnames(qs2) <- c("min", "q14", "median", "q34", "max")
is.nas2 <- df %>% is.na %>% as.data.frame %>% summarise_all(sum) %>% t()

stat2.1 <- data.frame(
  rownames=rownames(classes2),
  classes=classes2,
  is.nas=is.nas2,
  is.not.nas=nrow(df)-is.nas2
)

stat2.2 <- data.frame(
  rownames=rownames(means2),
  means=means2,
  vars=vars2,
  sds=sds2,
  qs2
)

stat2 <- full_join(stat2.1, stat2.2, by="rownames")
write.csv(x=stat2, file="summary2.csv")
```
|	|rownames	|classes	|is.nas	|is.not.nas	|means	|vars	|sds	|min	|q14	|median	|q34	|max|
|:--|:--|:--|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|
|1	|TrackID	|integer	|0	|33581	|16791	|93976428.5	|9694.144031	|1	|8396	|16791	|25186	|33581|
|2	|Study.Site	|character	|0	|33581	|NA	|NA	|NA	|NA	|NA	|NA	|NA	|NA|
|3	|Animal.code	|character	|0	|33581	|NA	|NA	|NA	|NA	|NA	|NA	|NA	|NA|
|4	|NZDT_datetime	|character	|0	|33581	|NA	|NA	|NA	|NA	|NA	|NA	|NA	|NA|
|5	|Latitude	|numeric	|0	|33581	|-44.55962314	|1.766917292	|1.329254412	|-45.30308	|-45.2932	|-45.28644	|-45.27773	|-42.12284|
|6	|Longitude	|numeric	|0	|33581	|170.2712649	|3.213393283	|1.792594009	|169.26576	|169.28312	|169.29043	|169.30543	|173.56041|
|7	|HDOP	|numeric	|0	|33581	|3.071507698	|3.44993206	|1.857399273	|0.8	|1.9	|2.6	|3.6	|12.7|
|8	|Season	|character	|0	|33581	|NA	|NA	|NA	|NA	|NA	|NA	|NA	|NA|

Latitude の q34, maxもマイナスの値となりちゃんと南半球を示しているのが分かる。平均座標を確認すると、

```R:コード
m2 <- leaflet() %>% addTiles()
m2 <- m2 %>% 
  addCircles(lng = means2[3], lat = means2[2]) %>% 
  setView(lng = means2[3], lat = means2[2], zoom = 4)
m2
```
![R_map.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/32063de8-7048-61c5-99d3-4df2ee3e6db5.png)

となり、ニュージーランドで測定されたデータに修正されている。

```R:コード
g2.Animal.code <- ggplot(df, aes(x = Animal.code))
g2.Animal.code <- g2.Animal.code + geom_bar()
g2.Animal.code <- g2.Animal.code + theme(axis.text.x = element_text(angle = 90, hjust = 1))
plot(g2.Animal.code)
```
![R_Animal_code.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/74d28cb2-786f-4b5d-1418-915e9b0c0d8c.png)

```R:コード
g2.Season <- ggplot(df, aes(x = Season))
g2.Season <- g2.Season + geom_bar()
plot(g2.Season)
```
![R_Season.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1077734/7faa7e6e-d064-f353-5283-90c44a47e205.png)

表記ゆれについても修正できているのが分かる。

# おわりに
今回基礎分析したのはデータの中に明らかなミスがあるデータセットであった。分析につかうデータが必ずしもきれいな状態であるとは限らないので、まずはデータの基礎分析をきっちりこなせるよう注意してデータに触れていきたいと思う。