udata <- read.csv("dalbo_2021a_meta_data.csv")
udata <- udata[which(as.numeric(udata$period) == 8),]

treatment <- rep(0, nrow(udata))
basin <- rep(0, nrow(udata))

for(ii in 1:nrow(udata)){
	aSS <- udata[ii, "aSS"]
	aSH <- udata[ii, "aSH"]
	aHS <- udata[ii, "aHS"]
	aHH <- udata[ii, "aHH"]

	aSS2 <- aSS - aHH
	aSH2 <- aSH - aHH
	aHS2 <- aHS - aHH
	aHH2 <- 0
	aSH3 <- aSH2 / aSS2
	aHS3 <- aHS2 / aSS2

	lambda <- -aSH3
	Lambda <- 1 - aHS3

	basin[ii] <- Lambda / (lambda + Lambda)
	treat <- aSH3 + (100 * aHS3)
	treatment[ii] <- treat
}

udata <- cbind(udata, treatment, basin)

papers <- unique(udata[,"paper"])

ntreat <- 0
for (ii in 1:NROW(papers)){
	zdata <- udata[which(udata$paper == papers[ii]),]
	treat <- unique(zdata[,"treatment"])
	ntreat <- ntreat + NROW(treat)
}

vstag <- rep(0, ntreat)
bstag <- rep(0, ntreat)
ll <- 0

for (ii in 1:NROW(papers)){
	zdata <- udata[which(udata$paper == papers[ii]),]
	treat <- unique(zdata[,"treatment"])
	#print(papers[ii])
	#print(treat)
	for(jj in 1:NROW(treat)){
		zzdata <- zdata[which(zdata$treatment == treat[jj]),]
		maxp <- max(zzdata[,"period"])
		print(paste(papers[ii], treat[jj], maxp))
		tstag <- 0
		tot <- 0
		for(kk in 1:nrow(zzdata)){
			tot <- tot + 1
			tstag <- tstag + as.numeric(zzdata[kk, "stag"])
		}
		ll <- ll + 1
		rstag <- tstag/tot
		vstag[ll] <- rstag
		bstag[ll] <- as.numeric(zzdata[1, "basin"])
		print(bstag[ll])
		print(vstag[ll])
			
	}
}

ddata <- read.csv("dalbo_2021a_new_data.csv")
ddata <- ddata[which(as.numeric(ddata$period) == 8),]

#C4530 and C6520

data1 <- ddata[which(ddata$treatment == "C4530"),]
data2 <- ddata[which(ddata$treatment == "C6520"),]
s1 <- mean(data1$action)
s2 <- mean(data2$action)
b1 <- 0.6
b2 <- 0.385
bstag <- append(bstag, b1)
bstag <- append(bstag, b2)
vstag <- append(vstag, s1)
vstag <- append(vstag, s2)



#dev.new()
jpeg("stag.jpg")
plot(bstag, vstag, xlim = c(0,1), ylim = c(0,1), xlab = "Basin of Stag", ylab = "Probability of Stag")
segments(0, 1/6, 1/6, 1/6)
segments(1/6,1/6, 1/6, 5/6)
segments(1/6,5/6, 1, 5/6)
segments(0,0, 0.5, 0, lty = 2)
segments(0.5,0, 0.5, 1, lty = 2)
segments(0.5,1, 1, 1, lty = 2)

dev.off()




