# 2. TYÖARKKI
#
# Aineistona "mammals" Weisbergin kirjasta "Regression Analysis "(1985, s. 144-5) 
#ja se on myös Venablesin ja Ripleyn (1994, 97, 99) kirjastossa MASS.
#Weisberg, S. (1985) Applied Linear Regression. 2nd edition. Wiley,
#     pp. 144-5.
#
#     Selected from: Allison, T. and Cicchetti, D. V. (1976) Sleep in
#     mammals: ecological and constitutional correlates. Science 194,
#    732-734.
#
# Tutkitaan aivojen koon riippuvuutta ruumiinpainosta nisäkkäillä.
# Tarkastellaan regressiota ja grafiikkaa erityisesti mallin diagnostisia
# tarkasteluja varten. Kysytään esimerkiksi, noudattavatko residuaalit
# normaalijakaumaa.
#
# sink("A:\\TA2.txt")
# komennolla voidaan tulostus ohjata tiedostoon "A:\TA2.txt"
#mammals <- read.table("A:\\TMallit\\Datat\\mammals.txt", header = TRUE, sep = ",")
library(mass)
data(mammals)
attach(mammals)
#data()       # tulostaa data-luettelon
#
# Katso komennolla "help(data=mammals)"
help(mammals) # tietoja datasta
species <- row.names(mammals)
species    	# tulostaa aineiston sisältamien lajien nimet
mammals		# tulostaa aineiston		
x <- body; y <- brain
plot(x,y)
identify(x,y,species)	# etsi hiirellä esimerkiksi ihminen
			# oik.p. painikkeella poistutaan kuviosta
######
postscript("mammals.eps")
plot(x,y)
dev.off() # Laite pois päältä
######
plot(log(x),log(y))
identify(log(x),log(y),species)	# etsi jälleen hiirellä joitain lajeja
species.lm <- lm(y ~ x)
summary(species.lm)
#####
postscript("mamres.eps")
par(mfrow=c(2,2))	# 2 saraketta ja 2 riviä kuvapohjaan
plot(x,y)
abline(species.lm)	# piirtää suoran pisteparveen
r <- species.lm$residuals
f <- species.lm$fitted.values # tallennetaan seuraavia tarkasteluja varten
qqnorm(r)
qqline(r)
# Tarkistetaan, noudattavatko residuaalit normaalijakaumaa
# Ilmeisesti eivät noudata. Mistä syystä?
lx <- log(x)
ly <- log(y)
species.llm <- lm(ly ~ lx)
summary(species.llm)
plot(lx,ly)
abline(species.llm)
rl <- species.llm$residuals
fl <- species.llm$fitted.values
qqnorm(rl)
qqline(rl)	# Mitäs sanot tästä kuviosta?
dev.off() # Laite pois päältä
########
par(mfrow=c(2,1))
plot(f,r)
hist(r)
plot(fl,rl); hist(rl) # Lisää diagnostisia tarkistuksia
# Kumpi kahdesta regressiosta sopii paremmin tilanteeseen?
mam.mat <- cbind(x,y,lx,ly) # "sarakkeesitojalla" muodostetaan matriisi
cor(mam.mat)                # korrelaatiomatriisi
round(cor(mam.mat),3)       # pyöristys helpottaa lukemista
par(mfrow=c(1,1))          # yksi kuvio pohjaa kohti
pairs(mam.mat)
# Tehdään kuviosta kopio. Tallennetaan ps-tiedostona.
postscript("A:\\TMallit\\Rscripts\\mamcor.ps",height=4)
# Tiedostoon tallennettu 4 tuumaa korkea kuvio ps-koodina
# Kuvio voidaan nyt liittää myöhemmin johonkin dokumenttiin
pairs(mam.mat)  # siirtää korrelaatiograafin tiedostoon
dev.off()       # kääntää laitteen "pois päältä, joten
plot(x,y)      # näkyy näytöllä
#Esim. Afrikan norsu (33) vaikuttava havainto aritm. asteikolla
#Poistetaan se ja katsotaan, miten tulokset muuttuvat
summary(lm(y~x))
summary(lm(y[-33]~x[-33]))
#Kaksi suurinta havaintoa ovat yhdessä vaikuttavia. Poistetaan ne
#ja tarkastellaan tulosta
x[33];y[33];x[19];y[19]
summary(lm(y[y<4603]~x[x<2547]))
q()