P13 Fabric Analysis.R

From gnewarchaeology wiki
Jump to: navigation, search


  • fabric_analysis.R

/!\ TODO: Syntax highlight + Literate programming with Org-mode

### p13 fabric_analysis ###
### workspace associato: fabric_analysis.RData ###

# load library
library(MASS)
library(boot)
library(CircStats)
library(RFOC)
library("heR.Misc")

# DATA ENTRY & TRANSFORMATION

## read data.frame
points <- read.csv("/home/dncgst/project/pirro13/manipulation_and_analysis/data_frame/points_all.csv")
records <- read.csv("/home/dncgst/project/pirro13/manipulation_and_analysis/data_frame/records_all_var.csv")
taph <- read.csv("/home/dncgst/project/pirro13/manipulation_and_analysis/data_frame/taph.txt")
record_selec <- read.csv("/home/dncgst/project/pirro13/manipulation_and_analysis/data_frame/records_selec.csv")

## s_index (indice di sfericità=l_max/l_med)
index <- records$l_max/records$l_med
records <- cbind(records,index)
class(records)

## fabric frame -na
indexmin <- which(records$index < 1.5)
fabric <- records[-indexmin,]
fabric <- fabric [-which(is.na(fabric$index)),]
class(fabric)
summary(fabric)

#! bearing variable = o in degrees
#! plunge_direction variable = p direction in degrees (n,s,e,w,etc...)
#! plunge variable = p inclination in degrees

## vector bear from fabric.frame where us=,column number 9
bear_a <-as.numeric(fabric[which(fabric$us == "A"),10])
bear_b <-as.numeric(fabric[which(fabric$us == "B"),10])
bear_c <-as.numeric(fabric[which(fabric$us == "C"),10])
bear_d <-as.numeric(fabric[which(fabric$us == "D"),10])

## vector bone from fabric.frame where us=,column number 9
bear_a_osso <-as.numeric(fabric[which(fabric$us == "A" & fabric$descrizione == "OSSO"),10])
bear_a_dente <-as.numeric(fabric[which(fabric$us == "A" & fabric$descrizione == "DENTE"),10])
bones_a <- c(bear_a_osso,bear_a_dente)
bear_b_osso <-as.numeric(fabric[which(fabric$us == "B" & fabric$descrizione == "OSSO"),10])
bear_b_dente <-as.numeric(fabric[which(fabric$us == "B" & fabric$descrizione == "DENTE"),10])
bones_b <- c(bear_b_osso,bear_b_dente)
bear_c_osso <-as.numeric(fabric[which(fabric$us == "C" & fabric$descrizione == "OSSO"),10])
bear_c_dente <-as.numeric(fabric[which(fabric$us == "C" & fabric$descrizione == "DENTE"),10])
bones_c <- c(bear_c_osso,bear_c_dente)
bear_d_osso <-as.numeric(fabric[which(fabric$us == "D" & fabric$descrizione == "OSSO"),10])
bear_d_dente <-as.numeric(fabric[which(fabric$us == "D" & fabric$descrizione == "DENTE"),10])
bones_d <- c(bear_d_osso,bear_d_dente)

## vector of NA*bear
na_bear_a <- which (is.na(bear_a))
na_bear_b <- which (is.na(bear_b))
na_bear_c <- which (is.na(bear_c))
na_bear_d <- which (is.na(bear_d))

## vector of NA*bones
na_bones_a <- which (is.na(bones_a))
na_bones_b <- which (is.na(bones_b))
na_bones_c <- which (is.na(bones_c))
na_bones_d <- which (is.na(bones_d))

## converts degrees to radians
bones_a_rad <- rad(bones_a)
bones_b_rad <- rad(bones_b)
bones_c_rad <- rad(bones_c)
bones_d_rad <- rad(bones_d)

## circular summary*bear (radians)
circ.summary(bear_a_rad)
circ.summary(bear_b_rad[-na_bear_b])
circ.summary(bear_c_rad[-na_bear_c])
circ.summary(bear_d_rad)
#! n=sample size, mean.dir=sample mean direction, rho=sample mean resultant length

# TESTS

## fabric!

## kuiper's test
## Performs Kuiper’s one-sample test of uniformity on the circle
kuiper(bear_a_rad, alpha=0.01)
kuiper(bear_b_rad, alpha=0.01)
kuiper(bear_c_rad, alpha=0.01)
kuiper(bear_d_rad, alpha=0.01)
#! Reject Null Hypothesis → uniformity of radians on the circle

## rayleigh's test
## Rayleigh Test of Uniformity: General Unimodal Alternative
r.test(bear_a, degree=TRUE)
r.test(bear_b[-na_bear_b], degree=TRUE)
r.test(bear_c[-na_bear_c], degree=TRUE)
r.test(bear_d, degree=TRUE)
#! r.bar=mean resultant length
#! alternative hypothesis is a unimodal distribution with unknown mean direction and unknown mean resultant length.

## watson's test
## Performs a Watson’s goodness of fit test for the von Mises or circular uniform distribution
watson(bear_a_rad, alpha=0.01, dist='uniform')
watson(bear_b_rad[-na_bear_b], alpha=0.01, dist='uniform')
watson(bear_c_rad[-na_bear_c], alpha=0.01, dist='uniform')
watson(bear_d_rad, alpha=0.01, dist='uniform')
#! Reject Null Hypothesis → circular uniform distribution
watson(bear_a_rad, alpha=0.05, dist='vm')
watson(bear_b_rad[-na_bear_b], alpha=0.05, dist='vm')
watson(bear_c_rad[-na_bear_c], alpha=0.05, dist='vm')
watson(bear_d_rad, alpha=0.05, dist='vm')

## bones!

## kuiper's test
## Performs Kuiper’s one-sample test of uniformity on the circle
kuiper(bones_a_rad, alpha=0.01)
kuiper(bones_b_rad, alpha=0.01)
kuiper(bones_c_rad, alpha=0.01)
kuiper(bones_d_rad, alpha=0.01)
#! Reject Null Hypothesis → uniformity of radians on the circle

## rayleigh's test
## Rayleigh Test of Uniformity: General Unimodal Alternative
r.test(bones_a, degree=TRUE)
r.test(bones_b[-na_bones_b], degree=TRUE)
r.test(bones_c[-na_bones_c], degree=TRUE)
r.test(bones_d, degree=TRUE)
#! r.bar=mean resultant length
#! alternative hypothesis is a unimodal distribution with unknown mean direction and unknown mean resultant length.

## watson's test
## Performs a Watson’s goodness of fit test for the von Mises or circular uniform distribution
watson(bones_a_rad, alpha=0.01, dist='uniform')
watson(bones_b_rad[-na_bones_b], alpha=0.01, dist='uniform')
watson(bones_c_rad[-na_bones_c], alpha=0.01, dist='uniform')
watson(bones_d_rad, alpha=0.01, dist='uniform')
#! Reject Null Hypothesis → circular uniform distribution
watson(bones_a_rad, alpha=0.05, dist='vm')
watson(bones_b_rad[-na_bones_b], alpha=0.05, dist='vm')
watson(bones_c_rad[-na_bones_c], alpha=0.05, dist='vm')
watson(bones_d_rad, alpha=0.05, dist='vm')

# DIAGRAM

## fabric!
rose.diag (bear_a, bins=4, prop=1,pts=T,shrink=1.2)
title(main="us A", line=3)
rose.diag (bear_b[-na_bear_b], bins=4, prop=1,pts=T,shrink=1.2)
title(main="us B", line=3)
rose.diag (bear_c[-na_bear_c], bins=4, prop=1,pts=T,shrink=1.2)
title(main="us C", line=3)
rose.diag (bear_d, bins=4, prop=1,pts=T,shrink=1.2)
title(main="us D", line=3)

x11()
m <- matrix(1:4, ncol=2, byrow=TRUE)
layout(m)
rose(bear_a, bins=4, rscale=NULL, labels=T, rings=TRUE, col="grey")
title(main="us A", line=3)
rose(bear_b[-na_bear_b], bins=4, rscale=NULL, labels=T, rings=TRUE, col="grey")
title(main="us B", line=3)
rose(bear_c[-na_bear_c], bins=4, rscale=NULL, labels=T, rings=TRUE, col="grey")
title(main="us C", line=3)
rose(bear_d, bins=4, rscale=NULL, labels=T, rings=TRUE, col="grey")
title(main="us D", line=3)

## bone!
x11()
m <- matrix(1:4, ncol=2, byrow=TRUE)
layout(m)
rose(bones_a, bins=32, rscale=NULL, labels=T, rings=TRUE, col="grey")
title(main="us A", line=3)
rose(bones_b[-na_bones_b], bins=32, rscale=NULL, labels=T, rings=TRUE, col="grey")
title(main="us B", line=3)
rose(bones_c[-na_bones_c], bins=32, rscale=NULL, labels=T, rings=TRUE, col="grey")
title(main="us C", line=3)
rose(bones_d, bins=32, rscale=NULL, labels=T, rings=TRUE, col="grey")
title(main="us D", line=3)

m <- matrix(1:4, ncol=2, byrow=TRUE)
layout(m)
rose.diag (bones_a, bins=4, prop=1,pts=T,shrink=1.2)
title(main="us A", line=3)
rose.diag (bones_b[-na_bones_b], bins=4, prop=1,pts=T,shrink=1.2)
title(main="us B", line=3)
rose.diag (bones_c[-na_bones_c], bins=4, prop=1,pts=T,shrink=1.2)
title(main="us C", line=3)
rose.diag (bones_d, bins=4, prop=1,pts=T,shrink=1.2)
title(main="us D", line=3)

## taphonomy!

abras <- which(taph$abrasione>=3)
select <- which(as.numeric (taph$ossid[abras]) >= 4)
#records_selec <- read.csv("/absolute/path/to/records_selec.csv")
#index <- records_selec$l_max/records_selec$l_med
#records_selec<- cbind(records_selec,index)
record_selec <- record_selec[select,]

select_a <- which(record_selec$us=="A")
na_a <- which (is.na(record_selec$bearing[select_a]))
select_b <- which(record_selec$us=="B")
na_b <- which (is.na(record_selec$bearing[select_b]))
select_c <- which(record_selec$us=="C")
na_c <- which (is.na(record_selec$bearing[select_c]))
select_d <- which(record_selec$us=="D")
na_d <- which (is.na(record_selec$bearing[select_d]))

#x11()
m <- matrix(1:2, ncol=2, byrow=TRUE)
layout(m)
rose.diag(record_selec$bearing[select_a[-na_a]], bins=4, prop=1,pts=T,shrink=1.2)
title(main="us A", line=3)
rose.diag (record_selec$bearing[select_b[-na_b]], bins=4, prop=1,pts=T,shrink=1.2)
title(main="us B", line=3)
rose.diag (record_selec$bearing[select_c], bins=4, prop=1,pts=T,shrink=1.2)
title(main="us C", line=3)
rose.diag (record_selec$bearing[select_d], bins=4, prop=1,pts=T,shrink=1.2)
title(main="us D", line=3)