## ----Household1, message=FALSE------------------------------------------------
mort <- agrodata::data_ibli("marsabit_mortality.rds")
mort[1:5, 1:7]


## ----loss1, message=FALSE-----------------------------------------------------
cause <- agrodata::data_ibli("marsabit_losses")
cause <- cause[,c("hhid", "sublocation", "cause","month", "year")]
colnames(cause)[4] <- "season"
s <- cause$season
cause$season <- "SRSD"
cause$season[s >=3 & s <=9] <-"LRLD"
cause <- na.omit(cause)
head(cause, n=2)



## ----loss2--------------------------------------------------------------------
mort <- merge(mort, cause, by=c("hhid","sublocation", "season", "year"))
mort <- mort[mort$cause == "Starvation/Drought", ]



## ----Household2---------------------------------------------------------------
mort <- mort[mort$type == "TLU", ]
amort <- aggregate(mort[,c("stock_beginning", "loss")], mort[ , c("year", "season", "sublocation")], sum, na.rm=TRUE)

head(amort, n=2)
amort$mortality_rate = 1 - ((amort$stock_beginning - amort$loss) / amort$stock_beginning )
amort$mortality_rate[amort$mortality_rate==0] <- NA


## ----hh2, message=FALSE-------------------------------------------------------
LR <- agrodata::data_ibli("marsabit_zLRLD")
SR <- agrodata::data_ibli("marsabit_zSRSD")

LR[1:3, 1:6]
lr <- reshape(LR, direction="long", varying = 2:ncol(LR), v.names="NDVI", timevar="year", times=2000:2015)
lr$id <- "LRLD"

sr <- reshape(SR, direction="long", varying = 2:ncol(LR), v.names="NDVI", timevar="year", times=2000:2015)
sr$id <- "SRSD"

ndvi <- rbind(lr, sr)
colnames(ndvi)[4] <- "season"
ndvi$year_season <- apply(ndvi[, c("year", "season")], 1, function(x) paste(x, collapse="_"))

head(ndvi, n=2)


## ----hh3----------------------------------------------------------------------
d <- merge(ndvi, amort, by=c("year", "season", "sublocation"))


## ----plots1-------------------------------------------------------------------
cols <- c("red", "blue")
seas <- c("SRSD", "LRLD")
pcols <- cols[match(d$season, seas)]

plot(mortality_rate~NDVI, data=d, xlab="zndvi", ylab="mortality", col=pcols, pch=20)

legend("topright", seas, col=cols, pch=20)



## ----r1, messsage=FALSE-------------------------------------------------------
dd <- na.omit(d[, c("NDVI", "mortality_rate")])
# add two "anchering points"
dd <- rbind(dd, c(-2.5, 0.9))
dd <- rbind(dd, c(-4, 1))
saveRDS(dd, "ndvi_mort.rds")

library(msir)
# Calculate and plot a 1.96 * SD prediction band
# that is a 95% prediction band
m <- loess.sd(dd)
plot(dd, cex=.5, col="blue", xlab="zNDVI", ylab="Mortality rate")
lines(m$x, m$y, col="red")
lines(m$x, m$upper, lty=2)
lines(m$x, m$lower, lty=2)


## ----ibli510------------------------------------------------------------------
x <- dd
x$NDVI = round(x$NDVI, 1)
a <- aggregate(x[, "mortality_rate", drop=FALSE], x[, "NDVI", drop=FALSE], mean)
ma <- loess.sd(a)
plot(a, cex=.5, col="blue", xlab="zNDVI", ylab="Mortality rate")
lines(ma$x, ma$y, col="red")
lines(ma$x, ma$upper, lty=2)
lines(ma$x, ma$lower, lty=2)


## ----ibli520------------------------------------------------------------------
ml <- loess(mortality_rate ~ NDVI, data=dd)
p <- predict(ml, d)
plot(d$mortality_rate, p, xlab="mortality rate", ylab="predicted")
abline(0,1)
x <- na.omit(cbind(d$mortality_rate, p))
cor(x[,1], x[,2])
d$predicted_mortality <- p
saveRDS(d, "pred_mort1.rds")


## ---- ibli530, fig.height=8---------------------------------------------------
ndvi$predicted_mortality <- predict(ml, ndvi)
par(mai=c(1,2,1,1))
boxplot(ndvi$predicted_mortality ~ ndvi$year_season, las=1, cex.axis=.5, horizontal=TRUE, xlab="zNDVI predicted Mortality rate", ylab="", cex=.25, col=rainbow(32))

saveRDS(ndvi, "pred_mort2.rds")