###################### Exercise sheet 3

#### Exercise 1 - Defining a binary indicator for the distance to a health hazard
## attach the asthma dataframe and load R packages 
library(splancs)
library(fields)
library(vcd)

load("asthma.RData.txt")
names(asthma)

## extract cases and controls and plot them
cases<-asthma[asthma$Asthma==1,]
controls<-asthma[asthma$Asthma==0,]
plot(controls$x,controls$y,xlab="xcoord",ylab="ycoord",main="Locations of cases and controls",col="green")
points(cases$x,cases$y,pch="+",col="red")

## extract a vector of distances to the chemical plant
dchem<-asthma$d2source2B
summary(dchem)

## define a binary indicator 
exp<-(dchem<=0.5)

## compute a contingency table 
ct<-table(exp,asthma$Asthma)
ct

## calculate the odds ratio by hand
or<-(ct[1,1]*ct[2,2])/(ct[1,2]*ct[2,1])
or

## print the odds ratio and the respective p-value
print(summary(oddsratio(ct)))
log(or)



###### Exercise 2 - Fitting a binomial GLM including distance to a health hazard to spatial case-control data
attach(asthma)

### fit a binomial GLM including the binary distance indicator
m1<-glm(Asthma~exp,family=binomial)
summary(m1)


### check the influence of different covariates who are suspected to increase the risk of asthma
m2<-glm(Asthma~exp+CAllergies,family=binomial)
summary(m2)

m3<-glm(Asthma~exp+CAllergies+Dust,family=binomial)
summary(m3)

m4<-glm(Asthma~CAllergies,family=binomial)
summary(m4)



####### Exercise 3 - Using a nonlinear function to model the distance from a health hazard
## WINDOWS USERS: Install package tribble2.zip using Packages --> Install from local zip file
### load library tribble 2
library(tribble2)

## check the distance from source 1 - coking works
source1<-tribbleNG(ccflag=Asthma,rho=sum(Asthma==1)/sum(Asthma==0),distances=d2source1B,alpha=1,beta=1)
source1

## check the distance from source 2 - chemical plant
source2<-tribbleNG(ccflag=Asthma,rho=sum(Asthma==1)/sum(Asthma==0),distances=d2source2B,alpha=0,beta=1)
source2


# include covariates
source1.HayFever<-tribbleNG(ccflag=Asthma,rho=source1$rho,distances=d2source1B,alpha=source1$alpha,beta=source1$beta,covars=HayFever,thetas=0)
source1.HayFever

source1.CAllergies<-tribbleNG(ccflag=Asthma,rho=source1$rho,distances=d2source1B,alpha=source1$alpha,beta=source1$beta,covars=CAllergies,thetas=0)
source1.CAllergies

source1.HayFever.CAllergies<-tribbleNG(ccflag=Asthma,rho=source1$rho,distances=d2source1B,alpha=source1$alpha,beta=source1$beta,covars=cbind(HayFever,CAllergies),thetas=c(0,0))
source1.HayFever.CAllergies


#### Exercise 4 - Estimating the ratio of two kernel estimates to get a risk surface
## load needed R library and dataset
library(splancs)
data(southlancs)

## extract larynx and lung cancer cases
larynx<-southlancs[southlancs[,3]==1,]
lung<-southlancs[southlancs[,3]==0,]

## extract locations and plot them
loclarynx<-cbind(larynx[,1],larynx[,2])
loclung<-cbind(lung[,1],lung[,2])
par(mfrow=c(1,2))
plot(loclarynx,xlab="eastings",ylab="northings",main="Locations of larynx cancer cases")
plot(loclung,xlab="eastings",ylab="northings",main="Locations of lung cancer cases")

## use a kernel smoother to estimate the continuous risk surface of both diseases 
k1<-kernel2d(loclung,southlancs.bdy,h0=2000)
k2<-kernel2d(loclarynx,southlancs.bdy,h0=2000)

## save the plots
pdf(file="lung.pdf")
filled.contour(k1,color=terrain.colors,main="Lung cancer")
dev.off()

pdf(file="larynx.pdf")
filled.contour(k2,color=terrain.colors, main="Larynx cancer")
dev.off()

## compute the ratio of the two kernels and save it 
kr<-kernrat(loclarynx,loclung,southlancs.bdy,h1=1500,h2=1500)

pdf(file="ratio.pdf")
filled.contour(kr,color=terrain.colors, main="Ratio")
dev.off()