# Análisis de la Covarianza.


# El archivo rotura_fibrilar.dat contiene el número de sesiones de sobrecarga excéntrica 
# recibidas por un pacientes con rotura fibrilar, así como el tiempo de recuperación de 
# la misma. Además los pacientes pueden haber recibido microgeneración endógena guiada 
# o no. Queremos investigar la influencia que tienen el número de sesiones de sobrecarga 
# excéntrica y la microgeneración endógena guiada sobre el tiempo de recuperación de la 
# rotura fibrilar.


dat<-read.table("rotura_fibrilar.dat", head=T)
attach(dat)

plot(w1,y, xlab="número de sesiones de sobrecarga excéntrica", ylab="tiempo de recuperación")

res<-lm(y~w1, data=dat)


plot(w1,y, type="n", xlab="número de sesiones de sobrecarga excéntrica", ylab="tiempo de recuperación")
text(w1,y,as.character(meg))

plot(w1,y, xlab="número de sesiones de sobrecarga excéntrica", ylab="tiempo de recuperación", col=c("red","blue")[meg], pch=20)

 
res1<-lm(y~w1+meg+w1:meg) # o equivalentemente res1<-lm(y~w1*meg)

model.matrix(res1)

summary(res1)


res2<-lm(y~w1+meg)
summary(res2)


step(res1)


abline(24.039, -0.572, col="red")
abline(24.039-4.817, -0.572, col="blue", lty=2)


layout(matrix(1:4,2))
plot(res2)



# El archivo catedral.dat contiene la altura (x) y la longitud (y) de la nave central
# de catedrales medievales inglesas (ambas medidas son en pies. 1 pie=0.3048m). Algunas
# de ellas son de estilo Románico (r) y otras de estilo Gótico (g). Algunas catedrales 
# tienen partes en ambos estilos y aparecen dos veces en los datos. Queremos investigar
# cómo la longitud está relacionada con la altura para los dos estilos.

 
datos<-read.table("catedral.dat",head=T)
attach(datos)
 
plot(x,y,type="n", xlab="Altura", ylab="Longitud")
text(x,y,as.character(estilo))
 
res<-lm(y~x+estilo+x:estilo) # o equivalentemente  res<-lm(y~x*estilo)

model.matrix(res)

summary(res)

res1<-lm(y~x+estilo)
summary(res1)

step(res)

abline(44.298,4.712)
abline(44.298+80.393,4.712,lty=2)

layout(matrix(1:4,2))
plot(res1)



# En el archivo iris tenemos las medidas en centímetros de las variables:
# 1. Longitud y anchura de los sepalos (Sepal.Length, Sepal.Width, respectivamente).
# 2. Longitud y anchura de los petalos (Petal.Length, Petal.Width, respectivamente).
# de 50 flores de cada una de 3 Especies de lirios (Setosa, Versicolor y Virgínica).
# (Datos publicados por E. Anderson(1935) y R.A. Fisher(1936)).

# Se pretende conocer si la anchura y la longitud de los sépalos están relacionadas
# y si esta relación depende de la especie de lirio de que se trate. A qué
# conclusiones podemos llegar a partir de los datos recogidos.

data(iris)
attach(iris)

plot(Sepal.Length,Sepal.Width)


res<-lm(Sepal.Width~Sepal.Length)
summary(res)


plot(Sepal.Length,Sepal.Width, col=c("red", "blue", "green3")[Species], pch=20)


res1<-lm(Sepal.Width~Sepal.Length*Species)
summary(res1)

beta<-coefficients(res1)

abline(a=beta[1],b=beta[2], col="red", lty=2) # Recta ajustada para setosa
abline(a=beta[1]+beta[3],b=beta[2]+beta[5], col="blue", lty=3) # Recta ajustada para versicolor
abline(a=beta[1]+beta[4],b=beta[2]+beta[6], col="green3", lty=4) # Recta ajustada para virgínica


layout(matrix(1:4,2))
plot(res1)



# Se seleccionan al azar cinco casas recientemente vendidas en
# tres distintas zonas residenciales (A, B y C) de cierta ciudad. El
# precio de venta (y) se compara con el valor catastral de la
# propiedad (x). Los datos, donde se dan el precio de venta y el
# valor catastral en miles de euros, se encuentran en el archivo
# casas.dat. Ajustar una ecuación de regresión lineal y determinar
# si las pendientes para las tres zonas residenciales son las
# mismas.

datos<-read.table("casas.dat",head=T)
attach(datos)

plot(x,y,type="n", xlab="valor catastral", ylab="precio de venta")
text(x,y,as.character(zona))

res<-lm(y~x*zona)

model.matrix(res)

summary(res)

abline(a=31.3724,b=.2325)
abline(a=31.3724-128.8064,b=.2325+2.4112, lty=2)
abline(a=31.3724+1.8886,b=.2325+.6786,lty=3)

layout(matrix(1:4,2))
plot(res)


a<-seq(min(x[zona=="A"]),max(x[zona=="A"]),len=100)
b<-seq(min(x[zona=="B"]),max(x[zona=="B"]),len=100)
c<-seq(min(x[zona=="C"]),max(x[zona=="C"]),len=100)
xx<-c(a,b,c)
z<-factor(rep(c("A","B","C"),c(100,100,100)))


con<-predict(res,data.frame(x=xx,zona=z),interval="confidence")
pre<-predict(res,data.frame(x=xx,zona=z),interval="prediction")

plot(x,y,type="n", xlab="valor catastral", ylab="precio de venta",ylim=c(min(pre[,2]),max(pre[,3])))
text(x,y,as.character(zona))
lines(xx,pre[,2],col=2)
lines(xx,pre[,3],col=2)

lines(xx,con[,3],col=2,lty=2)
lines(xx,con[,2],col=2,lty=2)

lines(xx,pre[,1],col=4)