### Motor de cohete 
fuerza<-c(2158.70,1678.15,2316,2061.30,2207.5,1708.3,1784.7,2575,2357.9,2256.7,2165.2,2399.55,1779.8,2336.75,1765.3,2053.5,2414.4,2200.5,2654.2,1753.7)
edad<-c(15.5,23.75,8,17,5.5,19,24,2.5,7.5,11,13,3.75,25,9.75,22,18,6,12.5,2,21.5)
datos<-data.frame(edad,fuerza)

plot(datos)

model<-lm(fuerza~edad, datos)
summary(model)
## 
## Call:
## lm(formula = fuerza ~ edad, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -215.98  -50.68   28.74   66.61  106.76 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2627.822     44.184   59.48  < 2e-16 ***
## edad         -37.154      2.889  -12.86 1.64e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 96.11 on 18 degrees of freedom
## Multiple R-squared:  0.9018, Adjusted R-squared:  0.8964 
## F-statistic: 165.4 on 1 and 18 DF,  p-value: 1.643e-10
plot(datos, xlab = "Edad del lote", ylab="Fuerza de enlace",col="blue4",  pch = 16, cex = .7)
abline(model, col="darkred")

par(mfrow=c(2,2))
plot(model)

#estimacion de sigma^2: SSE/n-2
sum(model$residuals**2)/18
## [1] 9236.381
#deviance==SSE
deviance(model)/18
## [1] 9236.381
#intervalos de confianza
confint(model)
##                  2.5 %    97.5 %
## (Intercept) 2534.99540 2720.6493
## edad         -43.22338  -31.0838
#ANOVA 
anova(model)
## Analysis of Variance Table
## 
## Response: fuerza
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## edad       1 1527483 1527483  165.38 1.643e-10 ***
## Residuals 18  166255    9236                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#intervalos de confianza
conf_interval <- predict(model,data.frame(edad), interval="confidence",
                         level = 0.95)
#intervalos de prediccion
pred_interval <- predict(model, data.frame(edad), interval="prediction",
                         level = 0.95)

par(mfrow=c(1,1))
plot(edad,fuerza, xlab = "Edad del lote", ylab="Fuerza de enlace",  pch = 16, cex = .7)
abline(model, col="darkred")
lines(edad[order(edad)], conf_interval[,2][order(edad)], col="blue", lty=2)
lines(edad[order(edad)], conf_interval[,3][order(edad)], col="blue", lty=2)
lines(edad[order(edad)], pred_interval[,2][order(edad)], col="orange", lty=2)
lines(edad[order(edad)], pred_interval[,3][order(edad)], col="orange", lty=2)

##Usando ggplot
library(ggplot2)

ggplot(datos,aes(x=edad,y=fuerza))+
  geom_point(color="black")+
  geom_smooth(method="lm",se=T,color="darkred")+
  ggtitle("Cohete")+
  theme(plot.title = element_text(hjust = 0.5))+
  labs(x="Edad", y="Fuerza", subtitle="Intervalos de predicción vs. Intervalos de confianza")+
  annotate("text", x=20,y=2400,label=paste("y=",round(summary(model)$coef[1,1],3),"+",round(summary(model)$coef[2,1],3),"x"), col="darkred")+
  geom_line(aes(x=edad[order(edad)], y=conf_interval[,2][order(edad)], colour="darkblue"))+
  geom_line(aes(x=edad[order(edad)], y=conf_interval[,3][order(edad)],colour="darkblue"))+
  geom_line(aes(x=edad[order(edad)], y=pred_interval[,2][order(edad)], colour="darkorange") )+
  geom_line(aes(x=edad[order(edad)], y=pred_interval[,3][order(edad)], colour="darkorange"))+
  scale_color_manual(name = "Intervalos", labels = c("Confianza", "Predicción"),values = c("darkblue", "darkorange")) 
## `geom_smooth()` using formula 'y ~ x'