# GLM1/ Yleistetyistä lineaarisista malleista (Venables & Ripley, 1999,1997,1994)

# Tehtävä 1.19
##############
library(MASS)
options(contrasts=c("contr.treatment", "contr.poly"))

data(anorexia)
help(anorexia)
# anorexia
 
ax.1 <- glm(Postwt ~ Prewt + Treat + offset(Prewt),
            family = gaussian, data = anorexia)       # Huom. offset
summary(ax.1)


ax.0 <- glm(Postwt ~ Prewt + Treat,
family = gaussian, data = anorexia)
summary(ax.0)

# Miten tulkitset tulokset offset-muuttujan vaikutuksen?
attach(anorexia)
PP <- Postwt-Prewt
ax.2 <- lm(PP ~ Prewt + Treat)
summary(ax.0)
plot(Treat,PP)
plot(Prewt,PP)
plot(Prewt,Postwt)

##############
# Tehtävä 1.27
##############
plot(Postwt ~ Prewt, data = anorexia, las = 1,
          xlab = "Paino ennen hoitoa",
          ylab = "Paino hoidon jälkeen",
          pch = as.integer(anorexia$Treat),
          col = as.integer(anorexia$Treat),
          main = "Hoidon vaikutus anoreksiaan")

par(mfrow = c(3,1))
PP <- Postwt-Prewt
plot(Treat,PP)
plot(Prewt,PP)
plot(Prewt,Postwt)

plot(PP ~ Prewt, data = anorexia, las = 1,
          xlab = "Paino ennen hoitoa",
          ylab = "Painon muutos",
          pch = as.integer(anorexia$Treat),
          col = as.integer(anorexia$Treat),
          main = "Hoidon vaikutus anoreksiaan")

ax.lm1 <- lm(PP ~ Prewt, data = anorexia,
                subset = Treat == "Cont")
ax.lm2 <- lm(PP ~ Prewt, data = anorexia,
                subset = Treat == "CBT")
ax.lm3 <- lm(PP ~ Prewt, data = anorexia,
                subset = Treat == "FT")
ax.lm <- lm(PP ~ Prewt, data = anorexia)

Prewt <- seq(70, 95, 1)
lines(Prewt, predict(ax.lm1, list(Prewt = Prewt)), lty = 1, col = 1)
lines(Prewt, predict(ax.lm2, list(Prewt = Prewt)), lty = 2, col = 2)
lines(Prewt, predict(ax.lm3, list(Prewt = Prewt)), lty = 3, col = 3)

wt <- pretty(70,95)
axis(side=4, at=wt, srt=90)

par(mfrow = c(2,2))
plot(PP ~ Prewt, data = anorexia,subset = Treat == "Cont",xlim=c(70,95),ylim=c(-10,20),
	pch = as.integer(anorexia$Treat),
        col = as.integer(anorexia$Treat))

lines(Prewt, predict(ax.lm1, list(Prewt = Prewt)), lty = 2, 
	col = 2)
plot(PP ~ Prewt, data = anorexia,subset = Treat == "CBT",xlim=c(70,95),ylim=c(-10,20),
	pch = as.integer(anorexia$Treat),
        col = as.integer(anorexia$Treat))
lines(Prewt, predict(ax.lm2, list(Prewt = Prewt)), lty = 1, 
	col = 1)
plot(PP ~ Prewt, data = anorexia,subset = Treat == "FT",xlim=c(70,95),ylim=c(-10,20),
	pch = as.integer(anorexia$Treat),
        col = as.integer(anorexia$Treat))
lines(Prewt, predict(ax.lm3, list(Prewt = Prewt)), lty = 3, 
	col = 3)
# plot(Treat,PP)
plot(PP ~ Prewt, data = anorexia, las = 1,xlim=c(70,95),ylim=c(-10,20),
          xlab = "Paino ennen hoitoa",
          ylab = "Painon muutos",
          pch = as.integer(anorexia$Treat),
          col = as.integer(anorexia$Treat),
          main = "Hoidon vaikutus anoreksiaan")
lines(Prewt, predict(ax.lm, list(Prewt = Prewt)), lty = 4, 
	col = 4)
                

################