# produce artificial tsunami death data and production estimation
# based on class on 12/13 2019
# prepared by M.Okumura
#
# produce 50 region data by random numbers
# pop by gammma distribution
# waveold by normal distribution neglecting negative value
# wave by normal distibution with correlated with waveold
pop <- ceiling(rgamma(50, 10, rate=0.1)) 
waveold <-rnorm(50,10,5)
wave <- 0.5*waveold+0.5*rnorm(50,10,5)
waveold <- waveold -(waveold<0)*waveold
wave <- wave -(wave<0)*wave
#
plot(waveold,wave)
cor(waveold,wave)
stone <- 1/(1+exp(-0.1*waveold +1))>runif(50)
## real relationship
drate <- 1/(1+exp(-(-15+wave-2.5*stone)))
death <- rbinom(50,pop,drate)
plot(drate, death)
#
# simple linear regression ignoring data generation structure
rslt0 <-lm((death/(pop+1))~wave+stone)
summary(rslt0)
#
# generalized linear model of binomial distribution (logit link)
rslt1 <- glm(cbind(death,pop-death)~wave, family=binomial)
rslt2 <- glm(cbind(death,pop-death)~wave+stone, family=binomial)
summary(rslt1)
summary(rslt2)
#
# inverse weighting by etone existance
rslt3 <- glm(stone ~ waveold,family=binomial)
summary(rslt3)
ir=c("red","green")
bg=c(2,3)
plogit <- function(x)
plogit <- 1/(1+exp(1.24188-0.10695*(x)))
plot(waveold,stone, xlim=c(0,25),ylim=c(0,1),pch=as.numeric(stone),col=ir[stone+1]) 
curve(plogit, xlim=c(0,25), ylim=c(0,1),add=TRUE)
# produce weightings
pstone <-plogit(waveold)
cnt1 <- sum(1/pstone * stone)
cnt2 <- sum( 1/(1-pstone) * (1-stone))
wgt <- 1/pstone * stone /cnt1+ 1/(1-pstone)*(1-stone)/cnt2 
plot(waveold,wgt, xlim=c(0,30), pch=as.numeric(stone),col=ir[stone+1])
# generalized linear model with inverse weightings
rslt4 <- glm(cbind(death,pop-death)~wave+stone, family=binomial,weight = wgt)
summary(rslt4)