Este curso fue originalmente creado por el profesor Edgar Acuña de la Universidad de Puerto Rico.

Curso de Regresión

#*****************************************
#Laboratorio 1 de Regresion Aplicada:
# Regresion lineal simple en R
# Enero del 2010 
#*****************************************
#Leyendo el conjunto de datos mortalidad de la internet
muertes <- read.table("http://www.udec.cl/~jbustosm/muertes.csv",
header=TRUE, sep=";")
#Mostrando los datos
muertes
#Haciendo un plot de tasa de mortalidad versus porcentaje de inmunizacion
x=muertes$porc.inmuniz
y=muertes$tasa.mort
win.graph()
plot(x,y,xlab="porcentaje de inmunizacion", ylab="tasa de mortalidad")
title("Relacion de la tasa de mortalidad con el porcentaje de inmunizacion")
pais=muertes$nacion
text(x,y,labels=as.character(pais),cex=.65,col="blue",pos=3, srt=30)
#Haciendo el ajuste por minimos cuadrados
l1<-lsfit(x,y)
#Mostrando los resultados del ajuste minimo cuadratico
l1
#Imprimiendo un resultado mas corto del ajuste minimocuadratico
ls.print(l1)
#Trazando la linea de regresión sobre el plot de puntos
abline(l1)
alfa=l1$coeff[1]
beta=l1$coeff[2]
text(50,100,bquote(hat(y)==.(alfa)+.(beta)*x))
#Extrayendo las  observaciones anormales 11 y 12  y creando un 
#nuevo conjunto muertes1
muertes1<-muertes[-c(11,12),]
#Haciendo el ajuste por minimos cuadrados excluyendo las 
#observaciones anormales y ploteando la linea de regresion para el nuevo 
#conjunto de datos
x1=muertes1$porc.inmuniz
y1=muertes1$tasa.mort
l2<-lsfit(x1,y1)
win.graph()
plot(x1,y1,xlab="porcentaje de inmunizacion", ylab="tasa de mortalidad")
abline(l2)
alfa1=l2$coeff[1]
beta1=l2$coeff[2]
text(50,100,bquote(hat(y)==.(alfa1)+.(beta1)*x))
ls.print(l2)

# ********************************************
# Segundo laboratorio de R
# Inferencia en regresion lineal simple
# Febrero 2006
# *********************************************
# Leyendo el conjunto de datos  mortalidad
muertes <- read.table("http://www.udec.cl/~jbustosm/muertes.csv",
header=TRUE, sep=";")
# Calculo de la linea de regresion usando el comando lsfit
l1<-lsfit(muertes$porc.inmuniz,muertes$tasa.mort)
# Calculo de la linea de regresion usando el comando lm
l2<-lm(tasa.mort~porc.inmuniz,data=muertes)
l2
# resumiendo los resultados de la linea de regresion
summary(l2)
# Imprimiendo los coeficientes estimados
summary(l2)$coef
beta=summary(l2)$coef[2,1]
eebeta=summary(l2)$coef[2,2]
# Hallando el intervalo de confianza del 95% para la pendiente Beta
bint<-c(beta-qt(.975,18)*eebeta,beta+qt(.975,18)*eebeta)
bint
# analisis de varianza para regresion
anova(l2)
#Hallando el intervalo de confianza para la media y
# el intervalo de prediccion de Y para varios valores de X. 
porc.inmuniz<-c(79,80,89,107)
porc.inmuniz<-as.data.frame(porc.inmuniz)
# Intervalo de confianza para la media de Y
predict(l2,porc.inmuniz,se.fit=T,interval=c("confidence"),level=.99)
# Intervalo de prediccion para Y 
predict(l2,porc.inmuniz,se.fit=T,interval=c("prediction"),level=.95)
# Haciendo bandas de confianza y bandas de prediccion
porc.inmuniz=seq(0,100,.5)
porc.inmuniz=as.data.frame(porc.inmuniz)
limic=predict(l2,porc.inmuniz,interval=c("confidence"),level=.95)
limip=predict(l2,porc.inmuniz,interval=c("prediction"),level=.95)
limites=cbind(porc.inmuniz,limic[,2:3],limip[,2:3])
plot(limites[,1],limites[,2],xlab="porcentaje de inmunizacion",ylab="tasa de mortalidad",ylim=c(-100,300),type="l",col=2)
title("Bandas de confianza y prediccion")
points(limites[,1],limites[,3],type="l",col=2)
points(limites[,1],limites[,4],type="l",col=4)
points(limites[,1],limites[,5],type="l",col=4)
points(muertes$porc.inmuniz,muertes$tasa.mort)
bandas <- expression("intervalo media", "intervalo prediccion")
legend(-3, 1, bandas, lty=1, col=c(2,4), cex=.7)
abline(l2)





# **************************************************************************************
#Laboratorio 4: Analisis de residuales en regresion lineal simple
#Marzo 2008
# **************************************************************************************
# llamando al conjunto de datos mortalidad infantil
muertes <- read.table("http://www.udec.cl/~jbustosm/muertes.csv",
header=TRUE, sep=";")
# Ajustando el modelo de regresion
 l2<-lm(tasa.mort~porc.inmuniz,data=muertes)
# Listado de las componentes del objeto l2
 attributes(l2)
 # imprimiendo los residuales usuales
 l2$res
# imprimiendo un resumen de la regresion
 l3<-summary(l2)
# imprimiendo un listado de las componentes del objeto l3
 attributes(l3)
# Extrayendo la desviacion estandar estimada
 l3$sigma
# Calculando los residuales estandarizados
 l2$res/l3$sigma
# Calculando los residuales estudentizados (internamente)
 studresi<-rstandard(l2)
# Cotejando la normalidad de los residuales estudentizados
win.graph()
 
par(mfrow=c(1,3),oma=c(1,1,1,1))
 hist(studresi)
 boxplot(studresi,main="boxplot de residuales")
 qqnorm(studresi)
 qqline(studresi)
title("Cotejando normalidad de residuales",outer=TRUE)

win.graph()
par(mfrow=c(2,2),oma=c(1,1,1,1))
qqnorm(studresi, main="Plot de Normalidad")

qqline(studresi)
plot(studresi, main="Plot de residuales")
abline(h=0,col=2)
plot(l2$fitted,studresi)
title("residuales versus valores ajustados",cex=.5)
abline(h=0,col=2)

plot(muertes$porc.inmuniz,studresi, main="residuales versus predictora",cex=.5)
abline(h=0,col=2)
title("Analisis de residuales",outer=TRUE)
# Construyendo la funcion para calcular el estadistico Durbin-Watson
dw<-function(e)
{
# eliminando el primer dato del vector original
e1<-e[-1] 
# eliminando el ultimo dato del vector original
e2<-e[-length(e)]
# creando las diferencias
diff<-e1-e2
dw<-sum(diff^2)/sum(e^2)
dw
}
#Leyendo el archivo de datos de mortalidad directamente de la internet
muertes <- read.table("http://www.udec.cl/~jbustosm/muertes.csv",
header=TRUE, sep=";")
attributes(muertes)
#Ajustando el modelo lineal
l1<-lm(tasa.mort~porc.inmuniz,data=muertes)
#Calculando el estadistico de Durbin Watson
dw1<-dw(l1$res)
cat("\n ","el estadistico Durbin Watson de la regresion lineal es=",dw1,"\n")

#***********************************************************
# Laboratorio  5. 
#Introduccion a Regresion lineal multiple
#************************************************************
#Leyendo el archivo de datos de millaje directamente de la internet
millaje <- read.table("	http://www.udec.cl/~jbustosm/millaje.csv",
header=TRUE, sep=";")
#millaje<-read.table(file="c:/millaje.txt",header=T)
attributes(millaje)
#Haciendo el plot matricial
win.graph()
pairs(millaje)
#haciendo el plot de mpg versus hp
win.graph()
plot(millaje$hp,millaje$mp)
#Extrayendo las variables del conjunto millaje
mpg<-millaje$mpg
hp<-millaje$hp
wt<-millaje$wt
sp<-millaje$sp
vol<-millaje$vol
#Ajustando el modelo lineal mpg vs hp
 l1<-lm(mpg~hp)
summary(l1)
#Ajustando  el modelo cuadratico de mpg versus hp
hp2<-hp^2
l2<-lm(mpg~hp+hp2)
summary(l2)
#Ajustando  el modelo hiperbolico de mpg versus hp
hp1<-1/hp
l3<-lm(mpg~hp1)
summary(l3)
#Ajustando  el modelo lineal de mpg versus las predictoras  hp y wt
l4<-lm(mpg~hp+wt)
summary(l4)
#Ajustando  el modelo lineal de mpg versus las predictoras  hp1 y wt
l5<-lm(mpg~hp1+wt)
summary(l5)
#Ajustando  el modelo lineal de mpg versus todas las predictoras
l6<-lm(mpg~vol+hp+sp+wt)
summary(l6)

#***********************************************************************
# Laboratorio 6.
# Calculo matricial en Regresion lineal multiple
#*********************************************************************
#Leyendo el archivo de datos de millaje directamente de la internet
millaje <- read.table("	http://www.udec.cl/~jbustosm/millaje.csv",
header=TRUE, sep=";")
#millaje<-read.table(file="c:/millaje.txt",header=T)
mdatos<-as.matrix(millaje)
dim(mdatos)
# Construyendo la  matriz X, cuya primera columna es de unos 
mx<-millaje[,c(2:5)]
col1<-rep(1,82)
mx<-cbind(col1,mx)
mx<-as.matrix(mx)
#calculo del número de observaciones y variables predictoras
n<-dim(mx)[1]
p<-dim(mx)[2]-1
# Calculando la transpuesta de mx
t(mx)
# Calculando la transpuesta de mx y mutiplicandola por mx
t(mx)%*%mx
#Calculo de la matriz sombrero
  h<-mx%*%solve(t(mx)%*%mx)%*%t(mx)
    #Calculo de los valores ajustados yhat=hat*y
  yhat<-h%*%mdatos[,1]
  #Calculo de los residuales
  resid<-mdatos[,1]-yhat
  #Calculo de la suma de cuadrados de los errores
  sse<-sum(resid^2)
  #calculo del cuadrado medio del error=mse
 df<-n-p-1
  mse<-sse/df
 #Calculo de la desviacion estandar estimada
  s<-mse^.5

  #***********************************************************************
  # Laboratorio  7:
  #   Inferencia en Regresion lineal multiple
  # 
  #************************************************************************
  #Leyendo el archivo de datos de millaje directamente de la internet
  millaje <- read.table("	http://www.udec.cl/~jbustosm/millaje.csv",
  header=TRUE, sep=";")
  #millaje<-read.table(file="c:/millaje.txt",header=T)
  l1<-lm(mpg~.,data=millaje)
  #Hallando los p-values para las pruebas de hipotesis
  summary(l1)
  anova(l1)
  #Modelo reducido para probar B(vol)=B(hP)=0
  l2<-lm(mpg~sp+wt, data=millaje)
  anova(l2)
  p=dim(millaje)[2]-1
  n=dim(millaje)[1]
  k=2
  #Suma de cuadrados de regresion del modelo completo
  a=sum(anova(l1)$Sum[-(p+1)])
  #Suma de cuadrados de regresion del modelo reducid
  b=sum(anova(l2)$Sum[-(k+1)])
  #Cuadrado Medio del error del modelo completo
  c=anova(l1)$Mean[p+1] 
  #Calculo del F parcial
  fp<-((a-b)/2)/c
  fp
  #Hallando el percentil de la F con alpha=.05
  qf(.95,k,n-p-1)
  # hallando el intervalo de confianza del 95% para el valor medio
  sp<-100
  wt<-20
  vol<-90
  hp<-50
  nuevo<-as.data.frame(cbind(sp,wt,vol,hp))
  nuevo
  predict.lm(l1,nuevo,se.fit=T,interval=c("confidence"),level=.95)
 #Hallando el ntervalo de prediccion del 99% para los mismos datos
predict.lm(l1,nuevo,se.fit=T,interval=c("prediction"),level=.99)

#***********************************************************************
# Laboratorio 9.  
#Diagnosticos de casos influenciales en Regresion 
#lineal multiple.
#Ejecutar primero la funcion acuinflu que aparece al final
#  
#*********************************************************************
#Leyendo el archivo de datos de millaje directamente de la internet
millaje <- read.table("	http://www.udec.cl/~jbustosm/millaje.csv",
header=TRUE, sep=";")
#millaje<-read.table(file="c:/millaje.txt",header=T)
# Hallando el modelo de regresion ajustado
l1<-lm(mpg~.,data=millaje)
#Hallando los residuales estudentizados internamente
rstint<-rstandard(l1)
rstint
#Hallando los residuales estudentizados externamente
rstext<-rstudent(l1)
rstext
#Determinado las observaciones que son outliers
which(abs(as.vector(rstint))>2)
#Ploteando los residuales y mostrando los outliers segun rsint
out1=rstint[abs(rstint)>2]
win.graph()
par(mfrow=c(2,2))
plot(rstint)
text(names(out1),out1,names(out1),adj=c(1,1))
#mostrando los outliers segun rstext
which(abs(as.vector(rstext))>2)
# Construyendo la  matriz X, cuya primera columna es de unos 
mx<-millaje[,c(2:5)]
col1<-rep(1,82)
mx<-cbind(col1,mx)
mx<-as.matrix(mx)
p<-dim(mx)[2]
ndatos<-dim(mx)[1]
# Hallando los valores leverages(diagonal de la matriz Hat) de cada observacion
leverages<-hat(mx)
#Determinando las observaciones que tienen un leverage alto
indlev=which(abs(leverages)>3*p/ndatos)
print(indlev)
#Ploteando los valores leverages
lev1=leverages[indlev]
#win.graph()
plot(leverages)
text(indlev,lev1,indlev,adj=c(1,1))
#calculando la distancia Cook de cada observacion
cookd<-cooks.distance(l1)
#win.graph()
plot(cookd)
#Determinando las observaciones influenciales segun distancia Cook
which(cookd>qf(0.50,p,ndatos-p))
#calculando los dffits
dfits<-dffits(l1)
#Determinando las observaciones influenciales segun dffits
which(abs(as.vector(dfits))>2*(p/ndatos)^.5)
dif1=dfits[abs(dfits)>2*(p/ndatos)^.5]
#win.graph()
plot(dfits)
text(names(dif1),dif1,names(dif1),adj=c(1,1))
#Calculado los dbfetas
dfb<-dfbetas(l1)
#Determinando las observaciones influenciales para cada uno de los coeficientes
for( i in 1:p)
{tempo<-which(abs(as.vector(dfb[,i]))>2*(ndatos)^(-.5))
cat("valores que influencian el coeficiente",i,"\n")
print(tempo)
}
#Calculando los covratios de cada observacion
cvrat<-covratio(l1)
#Determinando las observaciones influenciales segun los covratios
indcvr=which(abs(1-as.vector(cvrat))>3*p/ndatos)
print(indcvr)
cvr1=cvrat[indcvr]
win.graph()
plot(cvrat)
text(indcvr,cvr1,indcvr,adj=c(-1,-1),cex=.7)
#Determinando todas las observaciones influenciales con por lo menos uno
#de los criterios
#Nota: la funcion acuinflu es una version corregida de la funcion
#influence.measures disponible en R. 
tempo<-acuinflu(l1)
summary(tempo)
#****************************************************************************

acuinflu<-function (lm.obj) 
{#Esta funcion modifica la fucnion influence.measures disponible en R.
# Edgar Acuna, Marzo 2003
    is.influential <- function(infmat, n) {
        k <- ncol(infmat) - 4
        if (n <= k) 
            stop("Too few cases, n < k")
        absmat <- abs(infmat)
        result <- cbind(absmat[, 1:k] > 2/sqrt(n), absmat[, k + 1] > 
            3 * sqrt(k/n), abs(1 - infmat[, k + 2]) > (3 * 
            k)/n, infmat[, k + 3]>qf(0.5, k, n - k), 
            infmat[, k + 4] > (3 * k)/n)
        dimnames(result) <- dimnames(infmat)
        result
    }
    infl <- lm.influence(lm.obj)
    p <- lm.obj$rank
    e <- weighted.residuals(lm.obj)
    s <- sqrt(sum(e^2, na.rm = TRUE)/df.residual(lm.obj))
    xxi <- chol2inv(lm.obj$qr$qr, lm.obj$qr$rank)
    si <- infl$sigma
    h <- infl$hat
    dfbetas <- infl$coefficients/outer(infl$sigma, sqrt(diag(xxi)))
    vn <- variable.names(lm.obj)
    vn[vn == "(Intercept)"] <- "1_"
    colnames(dfbetas) <- paste("dfb", abbreviate(vn), sep = ".")
    dffits <- e * sqrt(h)/(si * (1 - h))
    cov.ratio <- (si/s)^(2 * p)/(1 - h)
    cooks.d <- ((e/(s * (1 - h)))^2 * h)/p
    dn <- dimnames(lm.obj$qr$qr)
    infmat <- cbind(dfbetas, dffit = dffits, cov.r = cov.ratio, 
        cook.d = cooks.d, hat = h)
    is.inf <- is.influential(infmat, sum(h > 0))
    ans <- list(infmat = infmat, is.inf = is.inf, call = lm.obj$call)
    class(ans) <- "infl"
    ans
}

#***********************************************************************
# Laboratorio 10.  Plots de residuales
#*********************************************************************
#Leyendo el archivo de datos de millaje directamente de la internet
millaje <- read.table("	http://www.udec.cl/~jbustosm/millaje.csv",
header=TRUE, sep=";")
#millaje<-read.table(file="c:/millaje.txt",header=T)
millaje
#Hallando la matriz de correlaciones
cor(millaje)
# Construyendo la  matriz X, cuya primera columna es de unos 
mx<-millaje[,c(2:5)]
col1<-rep(1,82)
mx<-cbind(col1,mx)
mx<-as.matrix(mx)
#Construyendo la matriz sombrero H
H<-mx%*%solve(t(mx)%*%mx)%*%t(mx)
#Hallando los residuales
n<-dim(mx)[1]
In=diag(0,n,n)
ehat<-(In-H)%*%millaje[,1]
#hallando los plots de residuales
win.graph()
par(mfrow=c(2,2))
plot(mx[,2],ehat,xlab="sp")
plot(mx[,3],ehat,xlab="wt")
plot(mx[,4],ehat,xlab="vol")
plot(mx[,5],ehat,xlab="hp")
#Hallando los residuales sin incluir la primera predictora(sp)
mx1<-mx[,-2]
H1<-mx1%*%solve(t(mx1)%*%mx1)%*%t(mx1)
ehat1<-(In-H1)%*%millaje[,1]
#Hallando los residuales de la primera predictora vs las otras
H1<-mx1%*%solve(t(mx1)%*%mx1)%*%t(mx1)
ehat1r<-(In-H1)%*%mx[,2]
#Hallando los residuales de incluir solo la variable volumen
mvol<-mx[,c(1,4)]
hvol<- mvol%*%solve(t(mvol)%*%mvol)%*%t(mvol)
ehvol<-(In-hvol)%*%millaje[,1]
#Hallando los residuales de peso (wt) versus volumen
ehwt.vol<-(In-hvol)%*%mx[,3]
#Haciendo el plot de regresion parcial de wt versus vol
win.graph()
plot(ehwt.vol,ehvol)
#Hallando los residuales de peso (hp) versus volumen
ehp.vol<-(In-hvol)%*%mx[,5]
#Haciendo el plot de regresion parcial de wt versus vol
win.graph()
plot(ehp.vol,ehvol)
# Cotejando la suposicion de normalidad
l1<-lm(mpg~wt,data=millaje)
summary(l1)
rstint<-rstandard(l1)
win.graph()
 par(mfrow=c(1,3))
 hist(rstint)
 boxplot(rstint)
 qqnorm(rstint)
 qqline(rstint)
#aplicando pruebas noparametricas
shapiro.test(rstint)
ks.test(rstint,"pnorm",0,1)
#Cotejando si la varianza es constante
win.graph()
plot(l1$fitted,rstint)
# Extra: Graficando los puntos y la curva ajustada encima
millaje$z<-1/millaje$wt
l2<-lm(mpg~z,data=millaje)
summary(l2)
win.graph()
#plot(millaje$z,millaje$mpg)
puntosx<-seq(15,60,length=100)
puntosy<-l2$coef[1]+l2$coef[2]/puntosx
plot(millaje$wt,millaje$mpg)
lines(puntosx,puntosy)


#***********************************************************************
# Laboratorio 11.  Transformacion para estabilizar la varianza
#  
#*********************************************************************
#Leyendo el archivo de datos de millaje directamente de la internet
millaje <- read.table("	http://www.udec.cl/~jbustosm/millaje.csv",
header=TRUE, sep=";")
#millaje<-read.table(file="c:/millaje.txt",header=T)
l1<-lm(mpg~.,data=millaje)
# Agrupando la variable de respuesta mpg en intervalos
nclases<-nclass.scott(millaje$mpg)
clases<-cut(millaje$mpg,nclases)
dmilla<-cbind(millaje$mpg,clases)
#Hallando la potencia de la transformacion
medias<-rep(0,nclases)
vars<-rep(0,nclases)
for(j in 1:nclases)
{medias[j]<-mean(dmilla[dmilla[,2]==j,])
vars[j]<-var(as.vector(dmilla[dmilla[,2]==j,]))
}
lmeans<-log(medias)
lvars<-log(vars)
lsfit(lmeans,lvars)
# El lsfit indica que la varianza es proporcional a la media al cuadrado
# una transformacion logaritmica en la variable de respuesta es recomendada 
mpglog<-log(millaje$mpg)
millaje1<-cbind(millaje,mpglog)
l2<-lm(mpglog~sp+wt+vol+hp,data=millaje1)
summary(l1)
summary(l2)
# Considerando que la varianza es proporcional a la media al cubo
# una transformacion  h(y)=y^-0.5 es realizada 
mpg05<-millaje$mpg^-0.5
millaje2<-cbind(millaje,mpg05)
l3<-lm(mpg05~sp+wt+vol+hp,data=millaje2)
summary(l3)
win.graph()
par(mfrow=c(1,3))
plot(l1$fitted,rstandard(l1),main="sin Transformacion",ylim=c(-3,3))
abline(h=0)
plot(l2$fitted,rstandard(l2),main="var~e(y)^2,w=log(y)",ylim=c(-3,3))
abline(h=0)
plot(l3$fitted,rstandard(l3),main="Var~E(y)^3, w=y^-1/2",ylim=c(-3,3))
abline(h=0)


#***********************************************************************
# Laboratorio 12.  Transformacion de Box y Tidwell
#           en un modelo de lineal multiple.
#  
#*********************************************************************
#Leyendo el archivo de datos de millaje directamente de la internet
millaje <- read.table("	http://www.udec.cl/~jbustosm/millaje.csv",
header=TRUE, sep=";")
#millaje<-read.table(file="c:/millaje.txt",header=T)
# Hallando el modelo de regresion ajustado
l1<-lm(mpg~.,data=millaje)
summary(l1)
betas<-l1$coeff
#Excluyendo el intercepto del vector de betas
betas1<-betas[-1]
#hallando las nuevas variables zetas
z1<-millaje$sp*log(millaje$sp)
z2<-millaje$wt*log(millaje$wt)
z3<-millaje$vol*log(millaje$vol)
z4<-millaje$hp*log(millaje$hp)
millaje1<-cbind(millaje,z1,z2,z3,z4)
l2<-lm(mpg~.,data=millaje1)
betas2<-l2$coeff
gammas<-betas2[c(6:9)]
#Hallando los alfas
alfas<-(gammas/betas1)+1
alfas
#Creando las nuevas variables
sp1<-millaje1$sp^alfas[1]
wt1<-millaje1$wt^alfas[2]
vol1<-millaje1$vol^alfas[3]
hp1<-millaje1$hp^alfas[4]
#regresion con todas las variables transformadas
l3<-lm(millaje1$mpg~sp1+wt1+vol1+hp1)
summary(l3)
#
#Haciendo nuevamente la transformacion desde el inicio pero sin usar vol
millaje2<-millaje[,-4]
l11<-lm(mpg~.,data=millaje2)
summary(l11)
betas11<-l11$coeff
#Excluyendo el intercepto del vector de betas1
betas12<-betas11[-1]
#hallando las nuevas variables zetas
z11<-millaje2$sp*log(millaje2$sp)
z21<-millaje2$wt*log(millaje2$wt)
z31<-millaje2$hp*log(millaje2$hp)
millaje2<-cbind(millaje2,z11,z21,z31)
l21<-lm(mpg~.,data=millaje2)
betas22<-l21$coeff
gammas1<-betas22[c(5:7)]
#Hallando los alfas1
alfas1<-(gammas1/betas12)+1
alfas1
#Creando las nuevas variables
sp11<-millaje2$sp^alfas1[1]
wt11<-millaje2$wt^alfas1[2]
hp11<-millaje2$hp^alfas1[3]
l5<-lm(millaje2$mpg~sp11+wt11+hp11)
summary(l5)
#Haciendo la regresion con solo las dos variables significativas
l6<-lm(millaje2$mpg~wt11+hp11)
summary(l6)


#***********************************************************************
# Laboratorio 13.  Transformacion de Box y Cox para remediar la 
#           falta de normalidad.
#  
#*********************************************************************
#Leyendo el archivo de datos de millaje directamente de la internet
millaje <- read.table("	http://www.udec.cl/~jbustosm/millaje.csv",
header=TRUE, sep=";")
# Hallando el modelo de regresion ajustado
l1<-lm(mpg~.,data=millaje)
summary(l1)
#Cargando la libreria MASS que contiene la funcion boxcox
library(MASS)
#Haciendo la transformacion
bc<-boxcox(l1,lambda=seq(-.6,.6,length=20),plotit=T)
#Transformando la respuesta
millaje1<-millaje
millaje1$mpg<-((millaje$mpg)^-0.22-1)/-0.22
l2<-lm(mpg~.,data=millaje1)
summary(l2)
# Viendo el efecto de la transformacion Box-Cox
rstint<-rstandard(l1)
rstint1<-rstandard(l2)
win.graph()
 par(mfrow=c(2,3))
 hist(rstint)
 boxplot(rstint)
title("Antes de la transformacion")
 qqnorm(rstint)
 qqline(rstint)
hist(rstint1)
 boxplot(rstint1)
title("Despues de la tranformacion")
 qqnorm(rstint1)
 qqline(rstint1)


#***********************************************************************
# Laboratorio 14.  Minimos cuadrados ponderados
#  
#*********************************************************************
#Leyendo el archivo de datos de millaje directamente de la internet
millaje <- read.table("	http://www.udec.cl/~jbustosm/millaje.csv",
header=TRUE, sep=";")
#millaje<-read.table(file="c:/millaje.txt",header=T)
#
#Extrayendo la primera y ultima observacion y consideramdo solo
#las columnas de variable de respuesta(1) y de wt(3).
#
millaje1<-millaje[-c(1,82),c(1,3)]
win.graph()
plot(millaje1[,2],millaje1[,1])
#Agrupando los valores de x
tabla<-table(millaje1[,2])
valsx<-dimnames(tabla)
valsx<-as.double(unlist(valsx))
nvalsx<-length(valsx)
#Calculando la varianza de cada grupo
varsg<-rep(0,nvalsx)
for(i in 1:nvalsx)
{grupox<-millaje1[millaje1[,2]==valsx[i],]
varsg[i]<-var(grupox[,1])
}
varsg
#Asignando los pesos a cada observacion (metodo I)
win.graph()
plot(varsg,valsx)
ndatos<-dim(millaje1)[1]
pesos<-rep(0,ndatos)
for(i in 1:ndatos)
{for(j in 1:nvalsx)
{if(millaje1[i,2]==valsx[j])
{pesos[i]<-1/varsg[j]}
}
}
pesos
#resultados de regresion sin ponderar
lw1<-lm(mpg~.,data=millaje1,weights=pesos)
summary(lw1)
#Asignando los pesos a cada observacion (metodo II)
valsx2<-valsx^2
s2<-lm(varsg~valsx+valsx2)
s2
valsx<-millaje1[,2]
valsx2<-millaje1[,2]^2
nuevo<-data.frame(valsx,valsx2)
pesos1<-predict.lm(s2,nuevo)
pesos1<-1/pesos1
pesos1
#resultados de regresion ponderada
lw2<-lm(mpg~.,data=millaje1,weights=pesos1)
summary(lw2)
#************************************************
#Aplicando regresion ponderada con la variable sp
#************************************************
nclases<-nclass.scott(millaje[,2])
clases<-cut(millaje[,2],nclases)
millaje2<-cbind(millaje[,1],clases)
#Calculando la varianza de cada grupo
varsg2<-rep(0,nclases)
for(i in 1:nclases)
{grupox<-millaje2[millaje2[,2]==i,]
varsg2[i]<-var(grupox[,1])
}
varsg2
#Asignando los pesos a cada observacion 
win.graph()
plot(varsg2,1:nclases)
nd<-dim(millaje)[1]
pesos2<-rep(0,nd)
for(i in 1:nd)
{for(j in 1:nclases)
{if(millaje2[i,2]==j)
{pesos2[i]<-1/varsg2[j]}
}
}
pesos2
#Resultados usando regresion ponderada
lw3<-lm(mpg~sp,data=millaje[,c(1,2)],weights=pesos2)
summary(lw3)
#Resultados sin usar regresion ponderada
lw4<-lm(mpg~sp,data=millaje)
summary(lw4)




#***********************************************************************
# Laboratorio 15. Regresion con variables cualitativas
#  
#*********************************************************************
#Leyendo el archivo de datos de millaje directamente de la internet
pesobebe <- read.table("http://www.udec.cl/~jbustosm/pesobebe.csv",
header=TRUE, sep=";")

#pesobebe<-read.table(file="c:/pesobebe.txt",header=T)
pbebe<-pesobebe[,10]
pmama<-pesobebe[,3]
fuma<-pesobebe[,5]
#regresion sonsiderando que la  relacion entre pesomama y pesobebe no
# es afectada por la variable fumar
l1<-lm(pbebe~pmama)
summary(l1)
#regresion considerando que fumar puede afectar la relacion pesomama-pesobebe
l2<-lm(pbebe~pmama+fuma) 
summary(l2)
#haciendo una regresion separada para cada grupo
pbebe0<-pbebe[fuma==0]
pmama0<-pmama[fuma==0]
win.graph()
plot(pmama0,pbebe0)
l3<-lm(pbebe0~pmama0)
summary(l3)
pbebe1<-pbebe[fuma==1]
pmama1<-pmama[fuma==1]
win.graph()
plot(pmama1,pbebe1)
l4<-lm(pbebe1~pmama1)
summary(l4)
#Probando la hipotesis de pendientes iguales
gl0<-summary(l3)$df[2]
gl1<-summary(l4)$df[2]
var0<-summary(l3)$sigma^2
var1<-summary(l4)$sigma^2
#Probando igualdad de varianzas
razon=var0/var1
if(razon<qf(.95,gl0,gl1))
{cat("se acepta igualdad")}
#Calculo de la varianza ponderada
varp<-(gl0*var0+gl1*var1)/(gl0+gl1)
varp
#**Calculo de las sumas de cuadrados sxx para cada grupo
ssx0<-sum(pmama0^2)-length(pmama0)*mean(pmama0)^2
ssx1<-sum(pmama1^2)-length(pmama1)*mean(pmama1)^2
#*Calculo de la prueba t
tcal<-(l3$coeff[2]-l4$coeff[2])/sqrt(varp*(1/ssx0+1/ssx1))
#Decision
if(abs(tcal)>qt(0.925,gl0+gl1))
{cat("Se rechaza hipotesis de igualdad")}
cat("Se acepta que las pendientes son iguales")




#***********************************************************************
# Laboratorio 16: Regresion logistica 
#  
#*********************************************************************
#Leyendo el archivo de datos de millaje directamente de la internet
pesobebe <- read.table("http://www.udec.cl/~jbustosm/pesobebe.csv",
header=TRUE, sep=";")
#pesobebe<-read.table(file="c:/bajopeso.txt",header=T)
# la primera columna de pesobebe contiene la variable bajopeso que se obtiene
#codificando la ultima variable peso. Si peso<2500 entonces bajopeso es 0 de lo
#contario vale 1
#Haciendo la regresion de bajopeso versus pesomama
l1=lm(bajopeso~pesomama,data=pesobebe)
#Ploteando bajopeso versus pesomama y la linea de regresion
win.graph()
plot(pesobebe$pesomama,pesobebe$bajopeso)
title("bajopeso versus pesomama")
abline(l1)
#Ploteando los residuales de la regresion
win.graph()
par(mfrow=c(1,2))
plot(l1$fitted,rstandard(l1))
abline(h=0)
plot(pesobebe$pesomama,rstandard(l1))
abline(h=0)
#Ploteando la curva logistica F(x)=1/(1+exp(-x)) en el rango (-10,10)
puntos<-seq(-10,10,length=100)
win.graph()
par(mfrow=c(1,1))
plot(puntos,1/(1+exp(-puntos)),type="l")
title("plot de la curva logistica en (-10,10)")
#Notar que hay algo de similtud con el plot anterior
#Agrupando la  variable pesomama en 5 intervalos para sacar estimaciones de 
#las probabilidades de bajopeso
ngrupos<-5 # el metodo de scott no conviene aqui
grupos<-cut(pesobebe[,3],ngrupos)
datos1<-cbind(pesobebe[,1],pesobebe[,3],grupos)
proby<-rep(0,ngrupos)
for(i in 1:ngrupos)
{datos2<-datos1[datos1[,3]==i,]
proby[i]<-sum(datos2[,1])/length(datos2[,1])
}
cat("las probabilidades estimadas en cada grupos son:",proby,"\n")
#Desafortunadamente los grupos no son homogeneos y no se prestan para
#aplicar regresion logistica en datos agrupados
#De ahora en adelante vamos a excluir la columna peso de pesobebe
pesobebe<-pesobebe[,-10]
# Haciendo la regresion logistica simple
logis1<-glm(bajopeso~pesomama,data=pesobebe,family=binomial)
summary(logis1)
#Haciendo la regresion logistica multiple
logis2<-glm(bajopeso~.,data=pesobebe,family=binomial)
summary(logis2)
#Haciendo otra vez la regresion logistica incluyendo solo las variables mas #significativas
logis3<-glm(bajopeso~pesomama+raza+fuma+hipertensio,data=pesobebe,family=binomial)
summary(logis3)
#Observando el valor de la devianza y del AIC el tercer modelo seria el mejor modelo
#sin embargo en lo que sigue vamos a considerar los resultados del segundo modelo
#********************************************************
#Prediciendo las clases con el segundo modelo
#Haciendo la clasificacion por el metodo mas simple comparando con p=0.5
#********************************************************
phat<-fitted.values(logis2)
nobs<-dim(pesobebe)[1]
clases<-rep(0,nobs)
for(i in 1:nobs)
{if(phat[i]>=0.5){clases[i]<-1}
}
errores<-sum(clases!=pesobebe[,1])
rate<-errores/nobs
cat("la tasa de mala clasificacion es=",rate,"\n")
#*************************************************
# Haciendo la clasificacion con el metodo mas complicado calculando la sensitividad y #especificidad
#********************************
p<-seq(.1,.9,length=9)
sensit<-rep(0,9)
especif<-rep(0,9)
for(j in 1:9)
{clases1<-rep(0,nobs)
for(i in 1:nobs)
{if(phat[i]>=p[j]){clases1[i]<-1}
}
pesobebe1<-cbind(pesobebe[,1],clases1)
sibajo<-pesobebe1[pesobebe1[,1]==1,]
nobajo<-pesobebe1[pesobebe1[,1]==0,]
sensit[j]<-mean(sibajo[,1]==sibajo[,2])
especif[j]<-mean(nobajo[,1]==nobajo[,2])
}
tabla<-cbind(p,sensit,especif)
cat("Sensitividad y especifidad para varios valores de p\n")
print(tabla)
#Haciendo el plot para hallar el p optimo
win.graph()
plot(p,sensit,type="l")
lines(p,especif)
text(p,sensit,labels=p)
title("Ploteando la sentividad y especificidad para varios p") 
# p=0.30 parece ser el optimo
#Ploteaando la curva ROC
win.graph()
plot(1-especif,sensit,type="l")
text(1-especif,sensit,labels=p)
title("La curva ROC")
#Notar que para p=.30 la curva esta mas cerca a la esquina superior izquierda 
#Clasificacion final
clasesf<-rep(0,nobs)
for(i in 1:nobs)
{if(phat[i]>=0.3){clasesf[i]<-1}
}
erroresf<-sum(clasesf!=pesobebe[,1])
ratef<-erroresf/nobs
cat("la tasa de mala clasificacion optima es=",ratef,"\n")

#***********************************************************************
# Laboratorio 17:  Seleccion de variables 
# a) Usando el metodo forward. Hace uso de la libreria leaps y de la 
#funcion regsubsets de dicha libreria
#  b) Seleccion de los mejores subconjuntos usando leaps y step con 
# el criterio AIC
# Usa la funcion selforw
# 
#*********************************************************************
#Leyendo el archivo de datos de grasa directamente de la internet
grasa <- read.table("http://www.udec.cl/~jbustosm/grasa.csv",
header=TRUE, sep=";")

#grasa<-read.table(file="c:/grasa.txt",header=F)
# El numero maximo de variables a entrar sera igual al numero de 
# predictoras del conjunto original
maxvar<-dim(grasa)[2]
#

#llamando a la libreria leaps
library(leaps)
#Aplicando el metodo forward
freg<-regsubsets(grasa~., data=grasa,method="forward",nvmax=maxvar)
#Mostrando la salida de todos los pasos con la estadisticas respectiva
selforw(grasa[,2:14],grasa[,1],.15)
#
# Aplicando el metodo de los mejores subconjuntos
#matrix  de predictoras
grasa.x<-grasa[,2:14]
#vector de respuesta
grasa.y<-grasa[,1]
#nombres de las variables predictoras
nombres<-colnames(grasa.x)
leaps(grasa.x,grasa.y,method="r2",nbest=2,names=nombres)
leaps(grasa.x,grasa.y,method="adjr2",nbest=2,names=nombres)
#Mejor modelo usando Cp de mallows
bcp<-leaps(grasa.x,grasa.y,method="Cp",nbest=1,names=nombres)
bcp$Cp
p<-2:maxvar
plot(p,bcp$Cp,type="l")
title("Gráfica del Cp de Mallows según el tamaño del modelo")
lines(2:maxvar,2:maxvar)
#Hallando el mejor subcjunto usando stepwise y el criterio AIC
l1<-lm(grasa~.,data=grasa)
step(l1,scope=~.,direction="backward")
#Hallando primero la regresion con la variable predictora mas correlacionada V7 

l2=lm(grasa~abdomen,data=grasa)
step(l2,scope=~.+edad+peso+altura+cuello+pecho+cadera+muslo+rodilla+tobillo+biceps+antebrazo+muneca,direction="forward")
######################################################################################################
backelim=function(x,y,alpha){
# Hace forward elimination using el paquete leaps 
# x es la matriz o data frame de variables independientes
# y es el vector de respuestas
require(leaps)
m<-ncol(x)                           # numero de variables independientes
n<-nrow(x)                           # tamano de la muestra
vm<-1:m
x=as.matrix(x)                    # en caso de que x sea data frame
pvmin<-rep(0,m)
regsubsets(x,y,method="backward")->out.x
pv.orig<-1-pf((out.x$rss[vm]-out.x$rss[vm+1])*(n-vm-1)/out.x$rss[vm+1],1,n-vm-1)
# sequential min of p-values down from full model
for (i in m:1){pvmin[i]<-min(pv.orig[m:i])}
sigma2=out.x$sserr/(n -1+ out.x$intercept - out.x$last)
cat("Eliminacion hacia atras",fill=T)
cat("",fill=T)
out=data.frame(p=c(vm,m+1),nvar=c(vm-1,m),rem.var=c(NA,colnames(x)[out.x$vorder-1]),
         pvmin=round(c(NA,pvmin),4),
          s=round(sqrt(out.x$rss/(n-c(vm,m+1))),3),
        r2=round(1-(out.x$rss/out.x$nullrss),3),
        r2adj=round(1-(out.x$rss/out.x$nullrss)*(n-1)/(n-c(vm,m+1)),3),
        Cp=round((out.x$rss/sigma2)-n+2*c(vm,m+1),3)
)
a=length(pvmin[pvmin>alpha])
cat("p=numero de coeficientes en el modelo",fill=T)
cat("nvar=p-1=numero de variables predictoras",fill=T)
cat("rem.var=la variable a ser removida, el modelo actual no incluye",fill=T)
cat("        esta variable",fill=T)
cat("pvmin=pvalue de la F parcial correspondiente a la variable menos importante en cada paso",fill=T)
cat("",fill=T)
out1=cbind(out[rev(((m+2)-a):(m+1)),1:4],out[rev(((m+1)-a):m),5:8])
print(out1)
}
#############################################################################
forwabic=function(x,y,alpha){
# Hace seleccion  forward usando el paquete leaps
# x es uma matriz o data frame de variables independentes
# y es el vector de respuestas
#alpha es el nivel de significacion para la f-parcial
# 
require(leaps)
m=ncol(x)            # numero de variables independientes
n=nrow(x)            # tamano de muestra
vm=1:m
x=as.matrix(x)      # si x es una data frame
pvmax<-rep(0,m)
out.x=regsubsets(x,y,nbest=2, method="forward")
pv.orig<-1-pf((out.x$rss[vm]-out.x$rss[vm+1])*(n-vm-1)/out.x$rss[vm+1],1,n-vm-1)
for (i in 1:m){pvmax[i]<-max(pv.orig[1:i])}  # sequential max of pvalues
cat("Seleccion Forward",fill=T)
cat("",fill=T)
out<-data.frame(p=c(vm,m+1),nvar=c(vm-1,m),add.var=c(NA,colnames(x)[out.x$vorder-1]),
         pvmax=round(c(NA,pvmax),4),
  aic=round(n * log(out.x$rss/n) +  2*c(vm,m+1),3),
   bic=round(n * log(out.x$rss/n) +  c(vm,m+1)*log(n),3),
   gcv=round(n*out.x$rss/((n-c(vm,m+1))^2),4)
)
a=length(pvmax[pvmax<alpha])
cat("p=numero de coeficientes en el modelo, p=1 es por el  intercepto",fill=T)
cat("nvar=p-1=numero de variables  predictoras",fill=T)
cat("add.var=la variable que ha sido anadida al modelo actual",fill=T)
cat("pvmax=p-value de F-parcial correspondiente a la variable mas importante en cada paso",fill=T)
cat("",fill=T)
print(out[2:(a+1),])
}
#################################################################
selforw=function(x,y,alpha){
# Hace seleccion  forward usando el paquete leaps
# x es uma matriz o data frame de variables independentes
# y es el vector de respuestas
#alpha es el nivel de significacion para la f-parcial
# 
require(leaps)
m=ncol(x)            # numero de variables independientes
n=nrow(x)            # tamano de muestra
vm=1:m
x=as.matrix(x)      # si x es una data frame
pvmax<-rep(0,m)
out.x=regsubsets(x,y,method="forward")
pv.orig<-1-pf((out.x$rss[vm]-out.x$rss[vm+1])*(n-vm-1)/out.x$rss[vm+1],1,n-vm-1)
for (i in 1:m){pvmax[i]<-max(pv.orig[1:i])}  # sequential max of pvalues
cat("Seleccion Forward",fill=T)
cat("",fill=T)
sigma2=out.x$sserr/(n -1+ out.x$intercept - out.x$last)
out<-data.frame(p=c(vm,m+1),nvar=c(vm-1,m),add.var=c(NA,colnames(x)[out.x$vorder-1]),
         pvmax=round(c(NA,pvmax),4),
        s=round(sqrt(out.x$rss/(n-c(vm,m+1))),3),
        r2=round(1-(out.x$rss/out.x$nullrss),3),
        r2adj=round(1-(out.x$rss/out.x$nullrss)*(n-1)/(n-c(vm,m+1)),3),
        Cp=round((out.x$rss/sigma2)-n+2*c(vm,m+1),3)
)
a=length(pvmax[pvmax<alpha])
cat("p=numero de coeficientes en el modelo, p=1 es por el  intercepto",fill=T)
cat("nvar=p-1=numero de variables  predictoras",fill=T)
cat("add.var=la variable que ha sido anadida al modelo actual",fill=T)
cat("pvmax=p-value de F-parcial correspondiente a la variable mas importante en cada paso",fill=T)
cat("",fill=T)
print(out[2:(a+1),])
}



#**************************************************************
#Laboratorio 18:
#Funciones para calcular Press y Validacion cruzada para Regresion 
#
#***************************************************************
PRESS=function (x) 
{#x es un objecto que sale de aplicar lm
    sum(resid(x)^2/(1 - lm.influence(x)$hat)^2)
}
cv10reg=function(data, folds=10,repet)
{#data: el conjunto de datos
#folds: numero de partes
#numero de repeticiones del experimento
#Suposicion: la primera columna contiene la variable de respuesta
n=dim(data)[1]
p=dim(data)[2]
nombres=colnames(data)
f1=as.formula(paste(nombres[1],".",sep="~"))
#print(f1)
EVC<-rep(0,repet)
for(i in 1:repet)
{
resid <- matrix(0, 1, folds)
azar <- data[rank(runif(n)),  ]
parti <- floor(n/folds)
for(j in 1:folds) {
cc <- ((j - 1) * parti + 1):(j * parti)
if(j == folds) {
cc <- ((j - 1) * parti + 1):n
}
datap <- azar[cc,  ]
# La muestra de prueba
datat <- azar[ - cc,  ]
#La muestra de entrenamiento
tempo = lm(f1, data = datat)
tempo1 <- predict(tempo, datap)
resid[j] <- sum((tempo1 - datap[, 1])^2)
}
EVC[i] <- sum(resid)/n
}
cat("Los estimados del error promedio de prediccion en cada prediccion son:\n")
print(EVC)
cat("La estimacion del error promedio de prediccion segun el numero de repeticiones dado sera\n")
EVC1<-mean(EVC)
EVC1
}

#************************************************
#Laboratorio 19: Diagnosticos de Multicolineralidad: VIF y
# Numero condicion y Regression Ridge
# Usa la libreria MASS y nuestra funcion acunaridge
#*************************************************
library(MASS)
# Ejemplo 1: Conjunto de datos millaje
millaje <- read.table("	http://www.udec.cl/~jbustosm/millaje.csv",
header=TRUE, sep=";")
# Hallando la matriz de correlacion de las predictoras
mcor<-cor(millaje[,2:5])
#Hallando los VIFs 
vif<-diag(solve(mcor))
cat("los VIF's son:\n")
vif
#recalculando los VIF's sin la variable hp
diag(solve(cor(millaje[,2:4])))
#Hallando el numero condicion
ev<-eigen(cor(millaje[,2:5]))
evals<-ev$values
evals
cond<-sqrt(evals[1]/ev$values[4])
cond

#***********************************************************
#Laboratorio 20: Regresion con Componentes Principales
# Requiere la libreria mva
#
#************************************************************
library(MVA)
# Ejemplo 1: Conjunto de datos millaje
millaje<-read.table(file="http://www.udec.cl/~jbustosm/millaje.csv",header=T, sep=";")
pcmilla<-prcomp(millaje[,2:5],retx=T,scale=T)
summary(pcmilla)
#Haciendo el screeplot para seleccionar los componentes a usar
screeplot(pcmilla)
#Haciendo regresion usando los componenentes principales
milla1<-as.data.frame(cbind(millaje$mpg,pcmilla$x))
colnames(milla1)[1]="mpg"
l1<-lm(mpg~PC1,data=milla1)
summary(l1)
l2<-lm(mpg~PC1+PC2,data=milla1)
# Notar que los coeficientes de regresion anteriores no varian
summary(l2)
l3<-lm(mpg~PC1+PC2+PC3,data=milla1)
summary(l3)
l4<-lm(mpg~PC1+PC2+PC3+PC4,data=milla1)
summary(l4)

#********************************************************
#Ejemplo 2:Conjunto de datos pollution 
pollution<-read.table("http://www.udec.cl/~jbustosm/pollution.csv",header=F, sep=";")
pcpollu<-prcomp(pollution[,1:15],retx=T,scale=T)
summary(pcpollu)
# Haciendo el screeplot para seleccionar las componentes a usar
screeplot(pcpollu,npcs=15,type="lines")
#Haciendo regresion usando los componenentes principales
pollu1<-as.data.frame(cbind(pollution$V16,pcpollu$x))
#colnames(pollu1)[1]="mort"
m1<-lm(V1~PC1,data=pollu1)
#Regresion stepwise usando los componentes principales y el criterio AIC
summary(m1)
step(m1,scope=~.+PC2+PC3+PC4+PC5+PC6+PC7+PC8+PC9+PC10+PC11+PC12+PC13+PC14+PC15,direction="forward")
m2<-lm(V1~.,data=pollu1)
step(m2,scope=~.,direction="backward")

# ****************************************************************
# Lababoratorio 21: Regresion usando arboles de decision
# Usa la libreria tree de Ripley
#****************************************************************
library(tree)
data(airquality)
airquality<-airquality[complete.cases(airquality), ] 
  arbol<-tree(Ozone~Solar.R+Temp,data=airquality)
  arbol
 summary(arbol)
  win.graph()
 plot(arbol, type="u")
  text(arbol)
mejorarbol<-prune.tree(arbol,best=5)
  mejorarbol
win.graph()
  plot(mejorarbol, type="u")
  text(mejorarbol)
gtemp<-seq(min(airquality$Temp),max(airquality$Temp),length=50)
 gradiation<-seq(min(airquality$Solar.R),max(airquality$Solar.R),length=50)
 grid<-cbind(gtemp,gradiation)
 grid1<-list(Solar.R=gradiation,Temp=gtemp)
 grid1<-expand.grid(grid1)
 estimado<-predict(arbol,grid1)
  grid2<-as.data.frame(grid1)
  matest<-matrix(estimado,50,50)
  persp(gradiation,gtemp,matest, theta=30, phi=45, xlab="radiation", ylab="temperature", zlab="ozone")