Primer ejemplo
Este será el primer ejemplo del modelo de regresión lineal simple usando R. Suponga que se cuenta con los siguientes datos sobre el cambio climático
################################################################
# Ejemplo1. Cambio climático #
################################################################
data<-read.table("climate-change-2016.txt", header=TRUE)
head(data)
## year Global_temp_anomaly CO2.ppm. DJIA
## 1 1959 0.0596 315.97 679.36
## 2 1960 0.0204 316.91 615.89
## 3 1961 0.0775 317.64 731.14
## 4 1962 0.0888 318.45 652.10
## 5 1963 0.1068 318.99 762.95
## 6 1964 -0.1495 319.62 874.13
tail(data)
## year Global_temp_anomaly CO2.ppm. DJIA
## 53 2011 0.5788 391.65 12217.60
## 54 2012 0.6240 393.85 13104.10
## 55 2013 0.6679 396.52 16576.66
## 56 2014 0.7408 398.65 17823.07
## 57 2015 0.8998 400.83 17425.03
## 58 2016 0.9363 404.21 19762.60
El objetivo del modelo es saber si con estos datos es posible predecir la temperatura global. Para ello primero veamos como se ven graficamente
plot(data$year,data$Global_temp_anomaly)
Y realicemos un rápido análisi exploratorio de las variables
#outliers
scatter.smooth(data$CO2.ppm.,data$Global_temp_anomaly)
boxplot(data$Global_temp_anomaly,sub=paste("Outliers: ", boxplot.stats(data$Global_temp_anomaly)$out))
boxplot.stats(data$Global_temp_anomaly)
## $stats
## [1] -0.14950 0.08880 0.29845 0.57830 0.93630
##
## $n
## [1] 58
##
## $conf
## [1] 0.1968963 0.4000037
##
## $out
## numeric(0)
boxplot(data$CO2.ppm.,sub=paste("Outlier rows: ", boxplot.stats(data$CO2.ppm.)$out))
#correlación
cor(data$Global_temp_anomaly,data$CO2.ppm.)
## [1] 0.947305
Para ajustar el modelo en R se ulitiza el siguiente código
#modelo
modelo<-lm(data$Global_temp_anomaly~data$CO2.ppm.)
modelo<-lm(Global_temp_anomaly~CO2.ppm., data)
print(modelo)
##
## Call:
## lm(formula = Global_temp_anomaly ~ CO2.ppm., data = data)
##
## Coefficients:
## (Intercept) CO2.ppm.
## -3.179330 0.009918
plot( data$CO2.ppm., data$Global_temp_anomaly)
abline(modelo)
Y es posible extraer la información que se considere importante
summary(modelo)
##
## Call:
## lm(formula = Global_temp_anomaly ~ CO2.ppm., data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19302 -0.06951 -0.01530 0.07302 0.17683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.1793300 0.1584344 -20.07 <2e-16 ***
## CO2.ppm. 0.0099179 0.0004482 22.13 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08848 on 56 degrees of freedom
## Multiple R-squared: 0.8974, Adjusted R-squared: 0.8956
## F-statistic: 489.7 on 1 and 56 DF, p-value: < 2.2e-16
#Estimación de la varianza
(summary(modelo)$sigma)**2
## [1] 0.007828923
#varianza de los coeficientes
summary(modelo)$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.179329991 0.1584343806 -20.06717 3.018799e-27
## CO2.ppm. 0.009917924 0.0004481656 22.13004 2.302800e-29
(summary(modelo)$coefficients[1,2])**2
## [1] 0.02510145
(summary(modelo)$coefficients[2,2])**2
## [1] 2.008524e-07
(summary(modelo)$coef[,2])^2
## (Intercept) CO2.ppm.
## 2.510145e-02 2.008524e-07
#qúe otra informacion tenemos?
names( summary(modelo) )
## [1] "call" "terms" "residuals" "coefficients"
## [5] "aliased" "sigma" "df" "r.squared"
## [9] "adj.r.squared" "fstatistic" "cov.unscaled"
hist(summary(modelo)$res)
Lo anterior usando el paquete base, pero también es posible usar paquetes con graficas mas avanzadas como ggplot
library(ggplot2)
library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
fit<-lm(Global_temp_anomaly~CO2.ppm., data)
ggplot(data,aes(x=Global_temp_anomaly,y=CO2.ppm.))+
geom_point(aes(Global_temp_anomaly,CO2.ppm.),color="darkblue")+
geom_smooth(method="lm",se=TRUE,color="darkred")+
scale_x_continuous(name="Temperatura global")+
scale_y_continuous(name="C02")+
ggtitle("Modelo para la temperatura global")+
annotate("text", x=0.15,y=390,label=paste("y=",round(summary(fit)$coef[1,1],3),"+",round(summary(fit)$coef[2,1],3),"x"))+
annotate("rect", xmin = 0.5, xmax =1, ymin = 340, ymax =350, fill="white", colour="red") +
annotate("text",x=0.75,y=345,label=paste("Var est= ",round((summary(modelo)$sigma)**2,4)))
## `geom_smooth()` using formula 'y ~ x'