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'