#########################33 #ejemplo 1.1 crea un mini modelo lineal y los grfica ########################### x=c(0.2,0.5,1,2,3) y=c(8,10,18,35,60) mod1=lm(y~x) summary(mod1) plot(x,y) lines(x,fitted(mod1)) points(x,fitted(mod1)) ##########################################3 #ejemplo 1.2 crea un modelo lineal abriendolo de unos datos de un FICHERO ############################################## dat=read.table(file="ejemplo1.2.txt") names(dat)=c("profundidad","oxigeno") mod=lm(oxigeno~profundidad,data=dat) summary(mod) x=dat$profundidad y=dat$oxigeno plot(dat) abline(mod) # ojo que hace esto en si points(dat$profundidad,fitted(mod)) anova(mod) ############################## #ejemplo 4.1 abriendo los datos de una base de datos previa ################################ library(faraway) # Carga el paquete Faraway, que contiene los ejemplos data(savings) # Carga la base de datos (del paquete Faraway) relativa a tasas de ahorro por países g=lm(sr~pop15+pop75+dpi+ddpi,savings) # Realiza la regresión lineal múltiple summary(g) # Resultados sobre los parámetros estimados, y el ajuste del modelo ############################################################################################################# #calculo detallado de los coeficientes y sus intervalos de confianza ############################################################################################################## beta=coef(g) # Devuelve los coeficientes del modelo x=model.matrix(g) # Devuelve la matriz del diseño n=nrow(x) # Devuelve el número de filas de una matriz o de un data.frame p=ncol(x) # Devuelve el número de columnas de una matriz o de un data.frame xtxi=solve(t(x)%*%x) # Cálculo de la matriz (X'X)^{-1} rss=deviance(g) # Devuelve la suma residual de cuadrados sigma2=rss/(n-p) # Estimación de la varianza del error errortipico=sqrt(sigma2*diag(xtxi)) # Cálculo de los errores típicos de los estimadores de los coeficientes niv=0.95 # Nivel de confianza t=qt(1-(1-niv)/2,n-p) # Cálculo del cuantil de la T de Student, necesario para los intervalos de confianza # para los coeficientes betainf=beta-t*errortipico # Extremos inferiores de los intervalos de confianza para los coeficientes betasup=beta+t*errortipico # Extremos superiores de los intervalos de confianza para los coeficientes covbeta=sigma2*xtxi # Matriz de varianzas-covarianzas de los estimadores de los coeficientes ############################################################################################################### #graficando la elipce de confianza de los paralos coeficientes ############################################################################################################### library(ellipse) # Carga el paquete ellipse para representar la elipse de confianza plot(ellipse(g,c(2,3)),type="l",xlim=c(-1,0)) # Representa la elipse de confianza para los coeficientes segundo y tercero points(beta[2],beta[3],pch=18) # Representa los coeficientes segundo y tercero, como un punto en el centro de la elipse points(0,0) # Representa el punto (0,0) para comprobar si pertenece a la elipse, # y en consecuencia si se acepta o se rechaza la hipótesis de que # ambos coeficientes valen cero. ############################################################################################################### # haciendo modelos de varias variables y comparando suma residual de cuadrado ############################################################################################################### g1575dpi=lm(sr~pop15+pop75+dpi,savings) # Modelo lineal con tres de las variables rss0=deviance(g1575dpi) # Suma residual de cuadrados del modelo con tres variables rss=deviance(g) # Suma residual de cuadrados del modelo con cuatro variables rss0-rss # Diferencia de sumas residuales de cuadrados anova(g1575dpi,g) # Test F del modelo simplificado, g1575dpi, #respecto del modelo más general, g. g75dpiddpi=lm(sr~pop75+dpi+ddpi,savings) # Modelo lineal con otras tres variables rss0=deviance(g75dpiddpi) # Suma residual de cuadrados del modelo con tres variables rss=deviance(g) # Suma residual de cuadrados del modelo con cuatro variables rss0-rss # Diferencia de sumas residuales de cuadrados anova(g75dpiddpi,g) # Test F del modelo simplificado, g75dpiddpi, respecto del # modelo más general, g. g15ddpi=lm(sr~pop15+ddpi,savings) # Modelo lineal con dos variables rss0=deviance(g15ddpi) # Suma residual de cuadrados del modelo con dos variables rss=deviance(g) # Suma residual de cuadrados del modelo con cuatro variables rss0-rss # Diferencia de sumas residuales de cuadrados f=((rss0-rss)/2)/(rss/(n-p)) # Estadístico F pvalue=1-pf(f,2,n-p) # Nivel crítico para el test F, calculado directamente anova(g15ddpi,g) # Test F del modelo simplificado, g15ddpi, respecto del modelo más general, g. g15=lm(sr~pop15,savings) anova(g15,g1575dpi) ############################################################################################################### # haciendo predicciones ############################################################################################################## predict(g) # Devuelve las estimaciones/predicciones para los individuos de la muestra predict(g,data.frame(pop15=30,pop75=2,dpi=1000,ddpi=5)) # Devuelve las estimaciones/predicciones para una nueva observación intest=predict(g,interval="confidence") # Devuelve las estimaciones /predicciones para los individuos de la muestra, junto # con intervalos de confianza para la estimación de la media condicionada intpred=predict(g,interval="prediction") # Devuelve las estimaciones /predicciones para los individuos de la muestra, junto # con intervalos de predicción predict(g,data.frame(pop15=c(30,40),pop75=c(3,1.5),dpi=c(1500,500),ddpi=c(2.5,4)),interval="prediction") # Devuelve # las estimaciones /predicciones para dos nuevas observaciones, junto # con intervalos de predicción cbind(intest,intpred[,2:3]) ########################################################################################################### ## 5.1. Diagnosis de observaciones atípicas o influyentes ########################################################################################################## library(faraway) # Carga el paquete Faraway, que contiene los ejemplos data(savings) # Carga la base de datos (del paquete Faraway) relativa a tasas de ahorro por países g=lm(sr~pop15+pop75+dpi+ddpi,savings) # Realiza la regresión lineal múltiple x=model.matrix(g) # Devuelve la matriz del diseño lev=hat(x) # Devuelve los leverages o apalancamientos (diagonal de la matriz hat) he=residuals(g) # Devuelve los residuos r=rstandard(g) # Devuelve los residuos estandarizados qqnorm(r) # Realiza el QQ-Plot de los residuos estandarizados qqline(r) # Traza la diagonal en el QQ-Plot windows() # Abre una nueva ventana para futuros gráficos, de modo que no desaparezca el que ya está creado t=rstudent(g) # Devuelve los residuos estudentizados d=cooks.distance(g) # Devuelve las distancias de Cook opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0)) # Abre una ventana, que va a recoger los cuatro gráficos plot(g) # Representa cuatro gráficos para la diagnosis de observaciones atípicas o influyentes par(opar) ################################## ## 5.2.1. Regresión polinómica ## ################################## library(faraway) # Carga el paquete Faraway, que contiene los ejemplos data(savings) # Carga la base de datos (del paquete Faraway) relativa a tasas de ahorro por países g1=lm(sr~ddpi,savings) # Regresión lineal simple sobre la variable ddpi attach(savings) # Carga las variables del "data frame" como vectores de memoria ddpi2=ddpi*ddpi # Construye una variable que contiene los cuadrados de ddpi g2=lm(sr~ddpi+ddpi2) # Regresión polinómica de grado dos (mediante regresión múltiple sobre ddpi y ddpi2) g2=lm(sr~ddpi+I(ddpi^2),savings) # Regresión polinómica de grado dos (igual que antes, pero sin necesidad de crear el vector ddpi2) g3=lm(sr~ddpi+I(ddpi^2)+I(ddpi^3),savings) # Regresión polinómica de grado tres g15eddpi2=lm(sr~pop15+I(exp(pop15))+ddpi+I(ddpi^2)) # Regresión múltiple, que mezcla funciones de dos variables pop15 y ddpi ########################### ## 5.2.2. Interacciones ## ########################### ca=c(750,800,850,900,950,1000,1050,1100,1150,1200) # Datos del ejemplo 5.2 d=c(280, 300,440,350,470,320, 380, 450, 500, 360) hueso=c(1.06,1.08,1.04,1.04,1.15,1.04,1.15,1.21,1.33,1.21) mod=lm(hueso~ca+d) # Regresión lineal múltiple, que incluye sólo los efectos principales, sin interacción summary(mod) cad=ca*d # Se crea un vector que contiene los productos de las dos variables explicativas continuas modi=lm(hueso~ca+d+cad) # Regresión lineal múltiple, con efectos principales e interacción modi=lm(hueso~ca+d+ca:d) # Se obtiene lo mismo, pero con una sintaxis que evita definir el vector con los productos modi=lm(hueso~ca*d) # Se obtiene lo mismo, pues ca*d incluye la interacción y los términos de orden inferior, dentro # de un planteamiento jerárquico summary(modi) ############################################################################################################### #ejemplo6.1 Ultimo sobre comparaciones simultaneas test de levene y FACTORES ############################################################################################################## rio=read.table(file="C:/Users/Anaderli/Desktop/Mod. Regresión/ejemplo6.1.txt",header=T) # Lectura de los datos del ejemplo mod=lm(oxigeno~grupo,rio) # Ejemplo (erróneo) de cómo saldría, si no se considera el grupo como factor xgrupo=factor(rio$grupo) # conversión del grupo en factor y=rio$oxigeno # El vector y recoge los valores de la variable respuesta modx=lm(y~xgrupo) # Modelo de análisis de la varianza #lm(oxigeno~grupo-1,rio) #abline(lm(oxigeno~grupo-1,rio)) anova(modx) # Tabla de análisis de la varianza y test F TukeyHSD(aov(y~xgrupo)) # Página 176 del libro de Faraway. Devuelve las comparaciones simultáneas abs_res=abs(residuals(modx)) # Devuelve los Valores absolutos de los residuos levene=lm(abs_res~xgrupo) # Realiza el test de Levene, mediante el test F de los valores absolutos de los residuos anova(levene) ##################################################################################################################### #ejemplo 7.1 library(MASS) # Carga la librería MASS (del libro de Venables y Ripley) data(whiteside) # Carga el ejemplo del consumo de gas, en función de la temperatura # y el aislamiento. Este ejemplo forma parte de la librería MASS attach(whiteside) # Permite usar las variables, sin poner delante el nombre del data frame ib=(Insul=="Before") # ib recoge los índices de los individuos que no tienen aislamiento ia=(Insul=="After") # ia recoge los índices de los individuos que tienen aislamiento plot(Temp[ib],Gas[ib]) # Representa un diagrama de dispersión de los individuos sin aislamiento points(Temp[ia],Gas[ia],pch=18) # Sobre el mismo diagrama, y con otros marcadores, # representa los individuos con aislamiento mod_0=lm(Gas~1) # Modelo que sólo contiene la constante, sin otra variable explicativa mod_t=lm(Gas~Temp) # Modelo de regresión del consumo de gas sobre la temperatura mod_i=lm(Gas~Insul) # Modelo de regresión del consumo de gas sobre el aislamiento mod_i_t=lm(Gas~Insul+Temp) # Modelo de regresión del consumo de gas sobre temperatura y aislamiento, # sin interacción mod_i_t_it=lm(Gas~Insul+Temp+Insul:Temp) # Modelo de regresión del consumo de gas sobre temperatura y aislamiento, # con interacción mod_i_t_it=lm(Gas~Insul*Temp) # Lo mismo, con otra sintaxis equivalente mod_t_ib=lm(Gas[ib]~Temp[ib]) # Modelo de regresión del consumo de gas sobre la temperatura, # para los individuos sin aislamiento abline(mod_t_ib) # Representa la recta de regresión del modelo anterior mod_t_ia=lm(Gas[ia]~Temp[ia]) # Modelo de regresión del consumo de gas sobre la temperatura, # para los individuos con aislamiento abline(mod_t_ia) # Representa la recta de regresión del modelo anterior # A partir de los coeficientes del modelo con interacción, se pueden calcular los coeficientes de las rectas # obtenidas separadamente para los grupos. # de la pifesora con interacciones tas.lm1<-lm(TAS~EDAD+SEXO+EDAD:SEXO, data=tension) #dentro del objeto, también se puede escribir la fórmula de forma abreviada: tas.lm1<-lm(TAS~EDAD*SEXO, data=tension) summary(tas.lm1)