library(arm)


# (a)
ceri=read.table(file.choose(),header=T)
str(ceri)

fit.1=glm(number~fuel+as.factor(strain),family=poisson,data=ceri)
display(fit.1)

# (b)
make.fake.poisson=function(x,fit)
{
	n=nrow(x)
	lambda=fitted(fit)
	x$number=rpois(n,lambda)
	return(x)
}
sum(ceri$number<=5)
#set.seed(85)
fake.count=replicate(1000,sum(make.fake.poisson(ceri,fit.1)$number<=5))
summary(fake.count)
histogram(fake.count,breaks=-.5:18.5)

# (c)
fit.2=glm(number~fuel+as.factor(strain),family=quasipoisson,data=ceri)
display(fit.2)
make.fake.quasipoisson=function(x,fit)
{
	n=nrow(x)
	mu=fitted(fit)
	alpha=mu/(1.2-1)
	x$number=rnbinom(n,mu=mu,size=alpha)
	return(x)
}
sum(ceri$number<=5)
#set.seed(85)
fake.count=replicate(1000,sum(make.fake.quasipoisson(ceri,fit.2)$number<=5))
summary(fake.count)
histogram(fake.count,breaks=-.5:18.5)



# (a)
hiv=read.table(file.choose(),header=T)
str(hiv)
attach(hiv)
xyplot(VL~GSS,data=hiv)


fit1=lm(VL~GSS)
coef(fit1)
se.coef(fit1)
sigma.hat(fit1)

# (b)
make.fake=function(x,fit)
{
	n=nrow(x)
	sigma=sigma.hat(fit)
	mu=fitted(fit)
	x$VL=rnorm(n,mu,sigma)
	return(x)
}
fit2=lm(VL~GSS,data=make.fake(hiv,fit1))
plot(fit2,which=1)
plot(fit1,which=1)

# (c)
fit3=lm(log(VL)~GSS)
coef(fit3)
se.coef(fit3)
sigma.hat(fit3)

make.fake=function(x,fit)
{
	n=nrow(x)
	sigma=sigma.hat(fit)
	mu=fitted(fit)
	x$VL=exp(rnorm(n,mu,sigma))
	return(x)
}
fit4=lm(log(VL)~GSS,data=make.fake(hiv,fit3))
plot(fit3,which=1)
plot(fit4,which=1)

# (d)
s=sim(fit3,1000)
pre=exp( s$beta %*% c(1,10) + rnorm(1000,0,s$gamma) )
quantile(pre,c(.025,.975))
exp( predict(fit3,data.frame(GSS=10)) )