.Rhistory

plot(model1,1,main="Model 1")
plot(model2,1,main="Model 2")
#create saturated model
model1.2 <- glm(cases~date+I(date^2)+I(date^3)+I(date^4),
data=aids,
family=poisson(link="log"))
#step until minimum AIC found
model1.3 <- step(model1.2)
model1final <- step(model1.3)
#do same for model 2
model2.2 <- glm(cases~date+I(date^2)+I(date^3)+I(date^4),
data=aids,
family=gaussian(link="log"))
model2.3 <- step(model2.2)
model2final <- step(model2.3)
#using anova to see if models are better
model1_anova <- anova(model1,model1final,test="Chisq")
model2_anova <- anova(model2,model2final,test="F")
print("Model 1 ANOVA: ")
model1_anova
print("Model 2 ANOVA: ")
model2_anova
par(mfrow=c(1,2))
plot(model1final,1,main="Model 1.2")
plot(model2final,1,main="Model 2.2")
#get preedicitons, se.fit will return standard errors
lambda_final <- predict(model1final,type="response",se.fit=TRUE)
aids$lambda_final <- lambda_final$fit
mu_final <- predict(model2final,type="response",se.fit=TRUE)
aids$mu_final <- mu_final$fit
#get confidence intervals, first lambda
lambda_final_se <- lambda_final$se.fit
lambda_final_z <- qnorm(1-0.05/2)
aids$lambda_final_lower <- aids$lambda_final - lambda_final_z*lambda_final_se
aids$lambda_final_upper <- aids$lambda_final + lambda_final_z*lambda_final_se
#same for mu
mu_final_se <- mu_final$se.fit
mu_final_z <- qnorm(1-0.05/2)
aids$mu_final_lower <- aids$mu_final - mu_final_z*mu_final_se
aids$mu_final_upper <- aids$mu_final + mu_final_z*mu_final_se
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))+
geom_line(aes(x=date,y=lambda_final,colour="Model 1 (Poisson)"),linewidth=1)+
geom_ribbon(aes(x=date,ymin=lambda_final_lower,ymax=lambda_final_upper,fill="Model 1 (Poisson)"),alpha=0.2)+
labs(x="Date",y="Cases",colour="Model",fill="Confidence Interval")+
scale_colour_manual(values=c("Model 1 (Poisson)"="blue"))+
scale_fill_manual(values=c("Model 1 (Poisson)"="blue"))+
guides(colour=guide_legend(order=1))
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))+
geom_line(aes(x=date,y=mu_final,colour="Model 2 (Normal)"),linewidth=1)+
geom_ribbon(aes(x=date,ymin=mu_final_lower,ymax=mu_final_upper,fill="Model 2 (Normal)"),alpha=0.2)+
labs(x="Date",y="Cases",colour="Model",fill="Confidence Interval")+
scale_colour_manual(values=c("Model 2 (Normal)"="red"))+
scale_fill_manual(values=c("Model 2 (Normal)"="red"))+
guides(colour=guide_legend(order=1))
library(MASS) #for glm.nb function
#make model
model3 <- glm.nb(cases~date+I(date^2)+I(date^3),
data=aids)
AIC(model3)
#get preedicitons, se.fit will return standard errors
binom <- predict(model3,type="response",se.fit=TRUE)
aids$binom <- binom$fit
#get confidence intervals, first lambda
binom_se <- binom$se.fit
binom_z <- qnorm(1-0.05/2)
aids$binom_lower <- aids$binom - binom_z*binom_se
aids$binom_upper <- aids$binom + binom_z*binom_se
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))+
geom_line(aes(x=date,y=binom,colour="Model 3 (Negative Binomial)"),linewidth=1)+
geom_ribbon(aes(x=date,ymin=binom_lower,ymax=binom_upper,fill="Model 3 (Negative Binomial)"),alpha=0.2)+
labs(x="Date",y="Cases",colour="Model",fill="Confidence Interval")+
scale_colour_manual(values=c("Model 3 (Negative Binomial)"="purple"))+
scale_fill_manual(values=c("Model 3 (Negative Binomial)"="purple"))+
guides(colour=guide_legend(order=1))
par(mfrow=c(1,2))
plot(model1,1,main="Model 1")
plot(model2,1,main="Model 2")
model2final$se.fit
mu_final$se.fit
BIC(model3)
AIC(model3)
BIC(model3)
BIC(model1)
BIC(model2_anova)
BIC(model2)
init_params <- c(100,0.1,-40)   #theta1, theta2, sigma
params <- nlm(f=mylike,
p=init_params,
x=nlmodel$x,
hessian=TRUE, #this will be used later to estimate errors
iterlim=10000, #increase iterations - more likely to converge
steptol=1e-10)
print("The estimates for theta1, theta2, and sigma are: ")
cat("theta1: ",params$estimate[1],"\n
theta2: ",params$estimate[2],"\n
sigma: ",params$estiamte[3])
init_params <- c(100,0.1,-40)   #theta1, theta2, sigma
params <- nlm(f=mylike,
p=init_params,
x=nlmodel$x,
hessian=TRUE, #this will be used later to estimate errors
iterlim=10000, #increase iterations - more likely to converge
steptol=1e-10)
cat("theta1: ",params$estimate[1],"\n
theta2: ",params$estimate[2],"\n
sigma: ",params$estiamte[3])
init_params <- c(100,0.1,-40)   #theta1, theta2, sigma
params <- nlm(f=mylike,
p=init_params,
x=nlmodel$x,
hessian=TRUE, #this will be used later to estimate errors
iterlim=10000, #increase iterations - more likely to converge
steptol=1e-10)
cat("theta1: ",params$estimate[1],"\n
theta2: ",params$estimate[2],"\n
sigma: ",params$estimate[3])
init_params <- c(100,0.1,-40)   #theta1, theta2, sigma
params <- nlm(f=mylike,
p=init_params,
x=nlmodel$x,
hessian=TRUE, #this will be used later to estimate errors
iterlim=10000, #increase iterations - more likely to converge
steptol=1e-10)
cat("theta1: ",params$estimate[1],
"theta2: ",params$estimate[2],
"sigma: ",params$estimate[3])
#find quartiles
cat(qnorm(0.995,0,1)
qnorm(0.005,0,1))
#find quartiles
cat(qnorm(0.995,0,1),
qnorm(0.005,0,1))
library(tidyverse)
library(ggplot2)
load("~/Rthur/StatisticalModelling/Coursework1/datasets_exercises.RData")
bad_model <- lm(y~x, data=nlmodel)
bad_data <- data.frame(x=seq(0,1,len=100))
bad_data$y <- predict(bad_model,newdata = bad_data)
ggplot(data=nlmodel,aes(x=x,y=y))+
geom_point()+
geom_line(data=bad_data,aes(x=x,y=y),linewidth=1,colour="red")
mylike <- function(parameters, x){
#set up values
y <- nlmodel$y
n <- length(y)
#unpack parameters
theta1 <- parameters[1]
theta2 <- parameters[2]
sigma <- parameters[3]
#find log likelihood
loglikelihood <- -50*log(2*pi)-50*log(sigma^2) - (1/(2*sigma^2))*sum((y - ((theta1 * x)/(theta2 + x)))^2)
return(-loglikelihood) #return negative (R minimises to optimise by default)
}
init_params <- c(100,0.1,-40)   #theta1, theta2, sigma
params <- nlm(f=mylike,
p=init_params,
x=nlmodel$x,
hessian=TRUE, #this will be used later to estimate errors
iterlim=10000, #increase iterations - more likely to converge
steptol=1e-10)
cat("theta1: ",params$estimate[1],
"theta2: ",params$estimate[2],
"sigma: ",params$estimate[3])
mean(nlmodel$y)
sd(nlmodel$y)
#get variances from the diagonal of the inverse hessian
variance_params <- diag(solve(params$hessian))
#square root to find stdev (standard error)
std_err_params <- sqrt(variance_params)
cat("Thus the standard error on theta1 is: ", std_err_params[1],
"and the standard error on theta2 is: ", std_err_params[2])
#find quartiles
cat(qnorm(0.995,0,1),
qnorm(0.005,0,1))
#unpack parameters
theta1 <- params$estimate[1]
theta2 <- params$estimate[2]
sigma <- params$estimate[3]
#theta1
theta1_confint <- c(theta1 - 2.576*std_err_params[1],
theta1 + 2.576*std_err_params[1])
#theta2
theta2_confint <- c(theta2 - 2.576*std_err_params[2],
theta2 + 2.576*std_err_params[2])
cat("The 99% confidence interval for theta1 is: ",theta1_confint,
"\nand the 99% confidence interval for theta2 is: ",theta2_confint)
zstat <- (params$estimate[2] - 0.08)/std_err_params[2]
zstat
z_extrema <- qnorm(c(0.05,0.95))
z_extrema
if (abs(zstat) > abs(z_extrema[1])){
cat("The z-statistic of ",zstat, "is outside of the significant region (",z_extrema,") in which we would reject the null hypothesis. Thus we accept the null hypothesis that theta2 != 0.08.")
} else {
cat("The z-statistic of ",zstat, "is inside the significant region(",z_extrema,"), meaning we reject the null hypothesis, and accept the new hypothesis, theta2 = 0.08.")
}
#create new data to predict over
xs <- seq(0,1,len=1000)
mus <- (theta1*xs)/(theta2+xs)
preds <- rnorm(n=length(xs),mean=mus,sd=sigma)
#95% prediction interval
#as model is normal we will use qnorm with our mu values and sigma
PI_upper <- qnorm(0.975,mean=mus,sd=sigma)
PI_lower <- qnorm(0.025,mean=mus,sd=sigma)
ggplot() +
geom_point(data = nlmodel, aes(x = x, y = y, colour = "data")) +
geom_line(aes(x = xs, y = mus, colour = "Mean Relationship"), linewidth = 1) +
geom_ribbon(aes(x = xs, y = mus, ymin = PI_lower, ymax = PI_upper, fill = "Prediction Interval"), alpha = 0.3) +
labs(x = "x", y = "y",
colour = "", fill = "") +
scale_colour_manual(values = c("data" = "black", "Mean Relationship" = "blue"),
labels = c("data" = "Data", "Mean Relationship" = "Mean Relationship")) +
scale_fill_manual(values = "blue", labels = c("Prediction Interval"))
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))
#fitting models
#both can be done with glm
#link="log" for both (we are using log of linear predictors lambda and mu)
model1 <- glm(cases~date,
data=aids,
family=poisson(link="log"))
model2 <- glm(cases~date,
data=aids,
family=gaussian(link="log"))
#get preedicitons, se.fit will return standard errors
lambda <- predict(model1,type="response",se.fit=TRUE)
aids$lambda <- lambda$fit
mu <- predict(model2,type="response",se.fit=TRUE)
aids$mu <- mu$fit
#get confidence intervals, first lambda
lambda_se <- lambda$se.fit
lambda_z <- qnorm(1-0.05/2)
aids$lambda_lower <- aids$lambda - lambda_z*lambda_se
aids$lambda_upper <- aids$lambda + lambda_z*lambda_se
#same for mu
mu_se <- mu$se.fit
mu_z <- qnorm(1-0.05/2)
aids$mu_lower <- aids$mu - mu_z*mu_se
aids$mu_upper <- aids$mu + mu_z*mu_se
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))+
geom_line(aes(x=date,y=lambda,colour="Model 1 (Poisson)"),linewidth=1)+
geom_line(aes(x=date,y=mu,colour="Model 2 (Normal)"),linewidth=1)+
geom_ribbon(aes(x=date,ymin=lambda_lower,ymax=lambda_upper,fill="Model 1 (Poisson)"),alpha=0.2)+
geom_ribbon(aes(x=date,ymin=mu_lower,ymax=mu_upper,fill="Model 2 (Normal)"),alpha=0.2)+
labs(x="Date",y="Cases",colour="Model",fill="Confidence Interval")+
scale_colour_manual(values=c("Model 1 (Poisson)"="blue","Model 2 (Normal)"="red"))+
scale_fill_manual(values=c("Model 1 (Poisson)"="blue","Model 2 (Normal)"="red"))+
guides(colour=guide_legend(order=1))
anova_result <- anova(model1,model2)
anova_result
#HMMM SHOULKD WE DO ANOVA OR DO AIC?
cat("AIC of Model 1: ",AIC(model1),"\n AIC of Model 2: ",AIC(model2))
par(mfrow=c(1,2))
plot(model1,1,main="Model 1")
plot(model2,1,main="Model 2")
#create saturated model
model1.2 <- glm(cases~date+I(date^2)+I(date^3)+I(date^4),
data=aids,
family=poisson(link="log"))
#step until minimum AIC found
model1.3 <- step(model1.2)
model1final <- step(model1.3)
#do same for model 2
model2.2 <- glm(cases~date+I(date^2)+I(date^3)+I(date^4),
data=aids,
family=gaussian(link="log"))
model2.3 <- step(model2.2)
model2final <- step(model2.3)
#using anova to see if models are better
model1_anova <- anova(model1,model1final,test="Chisq")
model2_anova <- anova(model2,model2final,test="F")
print("Model 1 ANOVA: ")
model1_anova
print("Model 2 ANOVA: ")
model2_anova
par(mfrow=c(1,2))
plot(model1final,1,main="Model 1.2")
plot(model2final,1,main="Model 2.2")
#get preedicitons, se.fit will return standard errors
lambda_final <- predict(model1final,type="response",se.fit=TRUE)
aids$lambda_final <- lambda_final$fit
mu_final <- predict(model2final,type="response",se.fit=TRUE)
aids$mu_final <- mu_final$fit
#get confidence intervals, first lambda
lambda_final_se <- lambda_final$se.fit
lambda_final_z <- qnorm(1-0.05/2)
aids$lambda_final_lower <- aids$lambda_final - lambda_final_z*lambda_final_se
aids$lambda_final_upper <- aids$lambda_final + lambda_final_z*lambda_final_se
#same for mu
mu_final_se <- mu_final$se.fit
mu_final_z <- qnorm(1-0.05/2)
aids$mu_final_lower <- aids$mu_final - mu_final_z*mu_final_se
aids$mu_final_upper <- aids$mu_final + mu_final_z*mu_final_se
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))+
geom_line(aes(x=date,y=lambda_final,colour="Model 1 (Poisson)"),linewidth=1)+
geom_ribbon(aes(x=date,ymin=lambda_final_lower,ymax=lambda_final_upper,fill="Model 1 (Poisson)"),alpha=0.2)+
labs(x="Date",y="Cases",colour="Model",fill="Confidence Interval")+
scale_colour_manual(values=c("Model 1 (Poisson)"="blue"))+
scale_fill_manual(values=c("Model 1 (Poisson)"="blue"))+
guides(colour=guide_legend(order=1))
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))+
geom_line(aes(x=date,y=mu_final,colour="Model 2 (Normal)"),linewidth=1)+
geom_ribbon(aes(x=date,ymin=mu_final_lower,ymax=mu_final_upper,fill="Model 2 (Normal)"),alpha=0.2)+
labs(x="Date",y="Cases",colour="Model",fill="Confidence Interval")+
scale_colour_manual(values=c("Model 2 (Normal)"="red"))+
scale_fill_manual(values=c("Model 2 (Normal)"="red"))+
guides(colour=guide_legend(order=1))
library(MASS) #for glm.nb function
#make model
model3 <- glm.nb(cases~date+I(date^2)+I(date^3),
data=aids)
AIC(model3)
#get preedicitons, se.fit will return standard errors
binom <- predict(model3,type="response",se.fit=TRUE)
aids$binom <- binom$fit
#get confidence intervals, first lambda
binom_se <- binom$se.fit
binom_z <- qnorm(1-0.05/2)
aids$binom_lower <- aids$binom - binom_z*binom_se
aids$binom_upper <- aids$binom + binom_z*binom_se
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))+
geom_line(aes(x=date,y=binom,colour="Model 3 (Negative Binomial)"),linewidth=1)+
geom_ribbon(aes(x=date,ymin=binom_lower,ymax=binom_upper,fill="Model 3 (Negative Binomial)"),alpha=0.2)+
labs(x="Date",y="Cases",colour="Model",fill="Confidence Interval")+
scale_colour_manual(values=c("Model 3 (Negative Binomial)"="purple"))+
scale_fill_manual(values=c("Model 3 (Negative Binomial)"="purple"))+
guides(colour=guide_legend(order=1))
#using anova to see if models are better
model1_anova <- anova(model1,model1final,test="Chisq")
model2_anova <- anova(model2,model2final,test="F")
model1_anova
model2_anova
#using anova to see if models are better
model1_anova <- anova(model1,model1final,test="Chisq")
model2_anova <- anova(model2,model2final,test="F")
model1_anova$`Pr(>Chi)`
model2_anova$`Pr(>F)`
#using anova to see if models are better
model1_anova <- anova(model1,model1final,test="Chisq")
model2_anova <- anova(model2,model2final,test="F")
model1_anova #$`Pr(>Chi)`
model2_anova #$`Pr(>F)`
#get preedicitons, se.fit will return standard errors
binom <- predict(model3,type="response",se.fit=TRUE)
aids$binom <- binom$fit
#get confidence intervals, first lambda
binom_se <- binom$se.fit
binom_z <- qnorm(1-0.05/2)
aids$binom_lower <- aids$binom - binom_z*binom_se
aids$binom_upper <- aids$binom + binom_z*binom_se
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))+
geom_line(aes(x=date,y=binom,colour="Model 3 (Negative Binomial)"))+
geom_ribbon(aes(x=date,ymin=binom_lower,ymax=binom_upper,fill="Model 3 (Negative Binomial)"),alpha=0.2)+
labs(x="Date",y="Cases",colour="Model",fill="Confidence Interval",title="Model 3 Predictions")+
scale_colour_manual(values=c("Model 3 (Negative Binomial)"="purple"))+
scale_fill_manual(values=c("Model 3 (Negative Binomial)"="purple"))+
guides(colour=guide_legend(order=1))
#get preedicitons, se.fit will return standard errors
binom <- predict(model3,type="response",se.fit=TRUE)
aids$binom <- binom$fit
#get confidence intervals, first lambda
binom_se <- binom$se.fit
binom_z <- qnorm(1-0.05/2)
aids$binom_lower <- aids$binom - binom_z*binom_se
aids$binom_upper <- aids$binom + binom_z*binom_se
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))+
geom_line(aes(x=date,y=binom,colour="Model 3 (Negative Binomial)"))+
geom_ribbon(aes(x=date,ymin=binom_lower,ymax=binom_upper,fill="Model 3 (Negative Binomial)"),alpha=0.2)+
labs(x="Date",y="Cases",colour="Model",fill="95% Confidence Interval",title="Model 3 Predictions")+
scale_colour_manual(values=c("Model 3 (Negative Binomial)"="purple"))+
scale_fill_manual(values=c("Model 3 (Negative Binomial)"="purple"))+
guides(colour=guide_legend(order=1))
#get preedicitons, se.fit will return standard errors
lambda_final <- predict(model1final,type="response",se.fit=TRUE)
aids$lambda_final <- lambda_final$fit
mu_final <- predict(model2final,type="response",se.fit=TRUE)
aids$mu_final <- mu_final$fit
#get confidence intervals, first lambda
lambda_final_se <- lambda_final$se.fit
lambda_final_z <- qnorm(1-0.05/2)
aids$lambda_final_lower <- aids$lambda_final - lambda_final_z*lambda_final_se
aids$lambda_final_upper <- aids$lambda_final + lambda_final_z*lambda_final_se
#same for mu
mu_final_se <- mu_final$se.fit
mu_final_z <- qnorm(1-0.05/2)
aids$mu_final_lower <- aids$mu_final - mu_final_z*mu_final_se
aids$mu_final_upper <- aids$mu_final + mu_final_z*mu_final_se
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))+
geom_line(aes(x=date,y=lambda_final,colour="Model 1 (Poisson)"))+
geom_ribbon(aes(x=date,ymin=lambda_final_lower,ymax=lambda_final_upper,fill="Model 1 (Poisson)"),alpha=0.2)+
labs(x="Date",y="Cases",colour="Model",fill="95% Confidence Interval",title="Model 1 Predictions")+
scale_colour_manual(values=c("Model 1 (Poisson)"="blue"))+
scale_fill_manual(values=c("Model 1 (Poisson)"="blue"))+
guides(colour=guide_legend(order=1))
ggplot(data=aids)+
geom_point(aes(x=date,y=cases))+
geom_line(aes(x=date,y=mu_final,colour="Model 2 (Normal)"))+
geom_ribbon(aes(x=date,ymin=mu_final_lower,ymax=mu_final_upper,fill="Model 2 (Normal)"),alpha=0.2)+
labs(x="Date",y="Cases",colour="Model",fill="95% Confidence Interval",title="Model 2 Predictions")+
scale_colour_manual(values=c("Model 2 (Normal)"="red"))+
scale_fill_manual(values=c("Model 2 (Normal)"="red"))+
guides(colour=guide_legend(order=1))
#import dependencies
library(ggplot2)
library(tidyverse)
library(fields)
load("datasets_project.RData")
df <- TBdata
head(df)
#-----Visualising Data------
plot.map(TBdata$TB[TBdata$Year==2024],n.levels=7,main="TB counts for 2014")
library(tidyverse)
library(fields)
install.packages("fields")
install.packages("maps")
install.packages("sp")
library(fields)
library(maps)
library(sp)
#get filepath
cwd <- getwd()
load("datasets_project.RData")
df <- TBdata
head(df)
#-----Visualising Data------
plot.map(TBdata$TB[TBdata$Year==2024],n.levels=7,main="TB counts for 2014")
#-----Visualising Data------
plot.map(TBdata$TB[TBdata$Year==2024],n.levels=7,main="TB counts for 2014")
#get filepath
cwd <- getwd()
load("datasets_project.RData")
df <- TBdata
head(df)
#-----Visualising Data------
plot.map(TBdata$TB[TBdata$Year==2024],n.levels=7,main="TB counts for 2014")
setwd("~/Rthur/StatisticalModelling/brazil_statistics")
#get filepath
cwd <- getwd()
load("datasets_project.RData")
df <- TBdata
head(df)
#-----Visualising Data------
plot.map(TBdata$TB[TBdata$Year==2024],n.levels=7,main="TB counts for 2014")
#-----Visualising Data------
plot.map(TBdata$TB[TBdata$Year==2024],n.levels=7,main="TB counts for 2014")
library(dplyr)
#-----Visualising Data------
plot.map(TBdata$TB[TBdata$Year==2024],n.levels=7,main="TB counts for 2014")
#-----Visualising Data------
plot.map(TBdata$TB[TBdata$Year==2014],n.levels=7,main="TB counts for 2014")
head(df)
model1 <- lm(TB~ns(Density,df=5)+ns(Poor_Sanitation,df=5)+Illiteracy,
data=TBdata)
library(spines)
install.packages("spines")
library(spines)
model1 <- lm(TB~ns(Density,df=5)+ns(Poor_Sanitation,df=5)+Illiteracy,
data=TBdata)
install.packages("splines2")
library(splines2)
model1 <- lm(TB~ns(Density,df=5)+ns(Poor_Sanitation,df=5)+Illiteracy,
data=TBdata)
library(splines)
model1 <- lm(TB~ns(Density,df=5)+ns(Poor_Sanitation,df=5)+Illiteracy,
data=TBdata)
plot.glm(model1)
plot.gam(model1)
sumary(model1)
summary(model1)
model1 <- lm(TB~ns(Density,df=5)+ns(Poor_Sanitation,df=5)+Illiteracy,
data=TBdata)
summary(model1)
plot.gam(model1)
library(mgcv)
plot.gam(model1)
plot.gam(model1)
?plot.gam
plot.gam(model,residuals=FALSE)
model1 <- gam(TB~Density+Poor_Sanitation)
summary(model1)
plot.gam(model,residuals=FALSE)
model1 <- gam(TB~Density+Poor_Sanitation)
model1 <- gam(TB~Density+Poor_Sanitation,
data=TBdata)
summary(model1)
plot.gam(model,residuals=FALSE)
model1 <- gam(formula=TB~Density+Poor_Sanitation,
data=TBdata)
summary(model1)
plot.gam(model,residuals=FALSE)
?gam
plot.gam(model,residuals=FALSE)
plot.gam(model1,residuals=FALSE)
plot.gam(model1)
summary(model1)
plot.gam(model1)
model1 <- gam(formula=TB~Density+Poor_Sanitation,
data=TBdata)
summary(model1)
plot.gam(model1)
#-----Visualising Data------
plot.map(TBdata$TB[TBdata$Year==2014],n.levels=7,main="TB counts for 2014")
head(TBdata)