setwd("D:/R Scripts/Linear Models/Data Files")
       library(MASS)

	par(mfrow = c(2., 2.), pty = "s", oma = c(0., 0., 4., 0.))
#	attach(cats)
      cats<-read.table("Cats.dat",header=T, col.names=c("Sex","Bwt","Hwt"))
      attach(cats)
	#-----------------------------------------------------
	#             LINEAR MODELS ON THE ORIGINAL SCALE
	#------------------------------------------------------
	plot(Bwt, Hwt, type = "n")
	text(Bwt, Hwt, c("+", "*")[Sex],col=c("red","blue"))
# Sex<-factor(Sex)
	cat("MODEL FOR BOTH SEXES\n\n")
	options(contrasts = c("contr.treatment", "contr.poly"))
	cats.fit <- lm(Hwt ~ Bwt + Sex + Bwt * Sex)
	print(summary(cats.fit, correlation = F))
	plot(Bwt, Hwt, type = "n")
	text(Bwt, Hwt, c("+", "*")[Sex])
	lines(Bwt[Sex == "F"], fitted(cats.fit)[Sex == "F"],col="red")
	lines(Bwt[Sex == "M"], fitted(cats.fit)[Sex == "M"],col="blue")
	cats.fit1 <- update(cats.fit,  ~ . - Bwt:Sex)
	print(anova(cats.fit1, cats.fit))
        stdres<-rstandard(cats.fit)
	par(mfrow = c(2., 2.))
	plot(fitted(cats.fit), stdres, type = "n")
	text(fitted(cats.fit), stdres, c("+", "*")[Sex])
	qqnorm(stdres)
	plot(fitted(cats.fit)[Sex == "F"], stdres[Sex == "F"], pch = "+")
	plot(fitted(cats.fit)[Sex == "M"], stdres[Sex == "M"], pch = "*")
	par(mfrow=c(2,1))
	plot(lm.influence(cats.fit)$hat, xlab = "index number", ylab = 
		"hat diagonal")
	abline(h = (2. * ncol(cats))/nrow(cats), lty = 2.)
	plot(cats.fit,which=4)
	#
	#-------------------------------------------------------
	#
	cat("\n\nSEPARATE MODEL FOR FEMALES\n")
	cats.f <- lm(Hwt ~ Bwt, subset = Sex == "F")
	print(summary(cats.f, correlation = F))
	cat("\n\nSEPARATE MODEL FOR MALES\n")
	cats.m <- lm(Hwt ~ Bwt, subset = Sex == "M")
	#--------Checking the equiality of variances------
	#
	#
	print(summary(cats.m, correlation = F))
	vf <- deviance(cats.f)/(nf <- cats.f$df.resid)
	vm <- deviance(cats.m)/(nm <- cats.m$df.resid)
	fstat <- c(Female = vf, Male = vm, F = (f <- vf/vm), Pr = pf(f, nf,
		nm))
	cat("\n\n\nF-test for equality of variances\n\n")
	print(fstat)
	par(mfrow = c(2, 2))
	boxcox(cats.fit)
	title("Box-Cox for Bwt")
	boxcox(lm(Hwt ~ log(Bwt) * Sex))
	title("Box-Cox for log(Bwt)")
	par(mfrow = c(2, 2))
	#-----------------------------------------------------
	#             LINEAR MODELS FOR LOG(RATIO)
	#------------------------------------------------------
	#
	lHwt <- log(Hwt)
	lBwt <- log(Bwt)
	lratio <- lHwt - lBwt
	cat("\n\nSEPARATE MODEL FOR FEMALES\n")
#	logcats.f <- lm(lratio ~ lBwt, subset = Sex == "F")
	logcats.f <- lm(lHwt ~ lBwt, subset = Sex == "F")
	#-------------------------------------------------------
	#
	print(summary(logcats.f, correlation = F))
	cat("\n\nSEPARATE MODEL FOR MALES\n")
#	logcats.m <- lm(lratio ~ lBwt, subset = Sex == "M")
	logcats.m <- lm(lHwt ~ lBwt, subset = Sex == "M")
	#
	print(summary(logcats.m, correlation = F))
	vf <- deviance(logcats.f)/(nf <- logcats.f$df.resid)
	vm <- deviance(logcats.m)/(nm <- logcats.m$df.resid)
	fstat <- c(Female = vf, Male = vm, F = (f <- vf/vm), Pr = pf(f, nf,
		nm))
	cat("\n\n\nF-test for equality of variances\n\n")
	print(fstat)
	cat("\n\nMODEL FOR LRATIO  (BOTH SEXES)\n")
#	logcats.fit <- lm(lratio ~ lBwt + Sex + lBwt:Sex, model = T)
	logcats.fit <- lm(lratio ~ lBwt + Sex + lBwt:Sex, model = T)
	#
	print(summary(logcats.fit, correlation = F))
	par(mfrow=c(2,2))
	plot(logcats.fit,which=c(1:2,4))
	plot(lm.influence(logcats.fit)$hat, xlab = "index number", ylab = 
		"hat diagonal")
	abline(h = (2. * ncol(cats))/nrow(cats), lty = 2.)
	par(mfrow = c(2., 2.))
	logcats.fin <- lm(lratio ~ 1.)
	#
	print(summary(logcats.fin, correlation = F))
# print(anova(logcats.fit,logcats.fin))
par(mfrow=c(2,2))
	plot(lBwt, lHwt, type = "n")
	text(lBwt, lHwt, c("+", "*")[Sex])
abline(a=logcats.fin$coefficients[1],b=1)
plot(lBwt,lratio,type="n")
text(lBwt, lratio, c("+", "*")[Sex])
abline(h=logcats.fin$coefficients[1])
	plot(Bwt, Hwt, type = "n")
	text(Bwt, Hwt, c("+", "*")[Sex])
abline(a=0,b=exp(logcats.fin$coefficients[1]))
plot(Bwt,Hwt/Bwt,type="n")
text(Bwt, Hwt/Bwt, c("+", "*")[Sex])
abline(h=exp(logcats.fin$coefficients[1]))
mtext("Final Model: lratio ~ 1",3,out=T)
par(mfrow=c(2,2))
	plot(logcats.fin,which=2)
plot(lBwt,logcats.fin$residuals)

#	alpha <- coefficients(logcats.fin)
	#	lratio.fit <- rep(alpha, length(lratio))
#	print(alpha)
#	lratio.fit <- fitted(logcats.fin)
#	resid <- lratio - lratio.fit
#	stdres <- resid/sqrt(var(resid))
#	plot(lBwt, lratio, type = "n")
#	text(lBwt, lratio, c("+", "*")[Sex])
#	lines(lBwt[Sex == "F"], lratio.fit[Sex == "F"])
#	lines(lBwt[Sex == "M"], lratio.fit[Sex == "M"])