ワークショップ「炭素年代測定値の分析と可視化」その1

7月19日に行われたワークショップ「炭素年代測定値の分析と可視化」の復習です。

当日資料はこちらのGitHubからダウンロードできます。
https://github.com/benmarwick/Analysing-Archaeological-Radiocabon-Ages-Using-R.git

単一年代の較正年代

厚真町厚幌2遺跡のデータを使用しました。
単一資料について較正年代を可視化します。

library(rcarbon)

x <- 
  calibrate(x = 4750,  # 1300BP
            errors = 30, # 誤差（標準偏差）
            calCurves = 'intcal13')  
# calCurvesは北半球、南半球、海産物などに応じて選択
# とりあえず"intcal13"を選択しておけば良い。

変数xはリスト形式のデータで、x$grids$1に較正年代とその密度が格納されています。

head(x$grids$"1")

      calBP       PrDens
44404  5597 1.431645e-05
44405  5596 2.038797e-05
44406  5595 2.874464e-05
44407  5594 4.010891e-05
44408  5593 5.550909e-05
44409  5592 7.618120e-05

plot(x,
     HPD=TRUE,
     credMass=0.95)

複数資料の較正年代

複数資料の較正年代も方法は上記と一緒です。
測定年代と誤差をベクトルで与えます。

xx <- 
  calibrate(x=c(4750,4910, 4950, 4960),
            errors=c(30,30, 30 ,30),
            calCurves='intcal13')

この場合、可視化手法は用意されていないので、summary()で概要を把握します。

summary(xx)

  DateID MedianBP OneSigma_BP_1 OneSigma_BP_2 TwoSigma_BP_1 TwoSigma_BP_2
1      1     5519  5580 to 5508  5486 to 5470  5585 to 5459  5373 to 5331
2      2     5632  5650 to 5605      NA to NA  5711 to 5687  5683 to 5674
3      3     5674  5716 to 5644  5625 to 5622  5733 to 5606      NA to NA
4      4     5684  5721 to 5652      NA to NA  5741 to 5609      NA to NA
  TwoSigma_BP_3
1      NA to NA
2  5663 to 5592
3      NA to NA
4      NA to NA

Bchronパッケージを利用した複数資料の分析と可視化

library(Bchron)

Bchronパッケージを使用して複数資料の年代の可視化を行います。

ages2 <- 
  BchronCalibrate(ages=c(4750,4910, 4950, 4960), # 測定年代
                  ageSds = c(30,30, 30 ,30),　# 誤差
                  calCurves= c('intcal13','intcal13','intcal13','intcal13')　# 較正方法
  )

plot(ages2)

ノンパラメトリックフェーズ推定

炭素年代から推定される、当該資料群採取遺跡（あるいは地域）の活動段階を推定する手法です。

# 引数はBchronCaribrate()と同じです
Dens <- 
  BchronDensity(ages=c(4750,4910, 4950, 4960), # 測定年代
                ageSds = c(30,30, 30 ,30),　# 誤差
                calCurves= c('intcal13','intcal13','intcal13','intcal13')　# 較正方法
  )

引数probの値を小さくするとより多くのphaseを検出します。phaseが多すぎると解釈が難しくなりそうです。

summary(Dens, prob = 0.95)

6つのフェーズが検出されました。

At probability level 0.95 there were 6 phases found.
Phase 1 
From 5317 to 5408 
Phase 2 
From 5411 to 5584 
Phase 3 
From 5596 to 5740 
Phase 4 
From 5757 to 5783 
Phase 5 
From 5787 to 5833 
Phase 6 
From 5837 to 5850 

この結果を較正年代とともに可視化します。

plot(Dens ,
     xlim = c(5300,5900) # x軸の範囲を指定
     )