Ejemplo hipotético de la relación entre la profundidad del agua (en centímetros) y el número de individuos de una rana. En este caso hipotético cuenta el número de ranas en 20 unidades de muestreo

¿Hay relación entre el número de ranas y la profundidad de las charcas?

datos del muestreo:

Profundidad_charco <- c(15,39,24,19,29,49,59,36,64,49,35,38,18,47,31,47,51,48,35,52)

No_ranas <- c(3,16,6,11,5,16,16,13,16,19,11,9,9,17,13,17,19,15,9,17)

# podemos visulalizar en forma de "tabla"
datos <- cbind(UM = rep(1:20), Profundidad_charco, No_ranas)

para graficar los datos:

par(mfrow= c(1,1))

plot(Profundidad_charco, No_ranas, pch = 16, col = "skyblue", cex = 2, bty = "l", xlim = c(0,69), ylim = c(0,30), las = 1, xlab = "Profundidad agua (cm)", ylab = "Número de ranas")

empleando un modelo lineal simple y se grafica la línea esperada

modelo <- lm(No_ranas ~ Profundidad_charco)
summary(modelo)

## 
## Call:
## lm(formula = No_ranas ~ Profundidad_charco)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.9729 -2.5913  0.4922  2.2025  3.8341 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         1.83270    1.96637   0.932    0.364    
## Profundidad_charco  0.28070    0.04741   5.921 1.33e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.842 on 18 degrees of freedom
## Multiple R-squared:  0.6607, Adjusted R-squared:  0.6419 
## F-statistic: 35.05 on 1 and 18 DF,  p-value: 1.327e-05

se grafica:

plot(Profundidad_charco, No_ranas, pch = 16, col = "skyblue", cex = 2, bty = "l", xlim = c(0,69), ylim = c(0,30), las = 1, xlab = "Profundidad agua (cm)", ylab = "Número de ranas")
abline(modelo, col= "red", lwd= 2)
text(15,28, "y = 1.83 + 0.28x", cex = 1.2)
text(15,25, "r2 = 0.66", cex = 1.2)

para entender las partes del modelo:

plot(Profundidad_charco, No_ranas, pch = 16, col = "skyblue", cex = 2, bty = "l", xlim = c(0,69), ylim = c(0,30), las = 1, xlab = "Profundidad agua (cm)", ylab = "Número de ranas")
abline(modelo, col= "red", lwd= 2)

text(60,28, "Media del modelo \n parte determinística")
arrows(60,26, 60, 19, length = 0.18, angle = 30)

segments(Profundidad_charco, fitted(lm(No_ranas ~ Profundidad_charco)), Profundidad_charco, No_ranas, lwd = 2, col = "black", lty = 3)

text(60,10, "Variación del modelo \n parte estocática")
arrows(51,10, 40, 9, length = 0.2, angle = 30)

Checando los supuestos:

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

Observando los datos extremos:

par(mfrow= c(1,1))
plot(Profundidad_charco, No_ranas, pch = 16, col = "skyblue", cex = 4, frame.plot = F, las = 1, xlab = "Profundidad agua (cm)", ylab = "Número de ranas")
abline(modelo, col= "red", lwd= 2)
text(Profundidad_charco, No_ranas, datos[,1])

Para guardar gráfico

jpeg(filename = "mi_primer_gráfico.jpg", width = 6000, height = 6000, units = "px", res =1200)

plot(Profundidad_charco, No_ranas, pch = 16, col = "skyblue", cex = 2, bty = "l", xlim = c(0,69), ylim = c(0,30), las = 1, xlab = "Profundidad agua (cm)", ylab = "Número de ranas")
abline(lm(No_ranas ~ Profundidad_charco), col= "red", lwd= 2)
text(20,28, "y = 1.83 + 0.28x", cex = 1.2)
text(20,25, "r2 = 0.66", cex = 1.2)

dev.off()

se pueden extraer los valores de los coeficientes y otros datos

#ejemplo el intercepto:

modelo$coefficients[1]

## (Intercept) 
##    1.832696

# la pendiente:
modelo$coefficients[2]

## Profundidad_charco 
##          0.2806957

# si se sugiere emplear la función "str" para ver detalles
str(modelo)

## List of 12
##  $ coefficients : Named num [1:2] 1.833 0.281
##   ..- attr(*, "names")= chr [1:2] "(Intercept)" "Profundidad_charco"
##  $ residuals    : Named num [1:20] -3.04 3.22 -2.57 3.83 -4.97 ...
##   ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:20] -57.47 -16.83 -1.38 5.32 -4.07 ...
##   ..- attr(*, "names")= chr [1:20] "(Intercept)" "Profundidad_charco" "" "" ...
##  $ rank         : int 2
##  $ fitted.values: Named num [1:20] 6.04 12.78 8.57 7.17 9.97 ...
##   ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
##  $ assign       : int [1:2] 0 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:20, 1:2] -4.472 0.224 0.224 0.224 0.224 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:20] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "(Intercept)" "Profundidad_charco"
##   .. ..- attr(*, "assign")= int [1:2] 0 1
##   ..$ qraux: num [1:2] 1.22 1.07
##   ..$ pivot: int [1:2] 1 2
##   ..$ tol  : num 1e-07
##   ..$ rank : int 2
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 18
##  $ xlevels      : Named list()
##  $ call         : language lm(formula = No_ranas ~ Profundidad_charco)
##  $ terms        :Classes 'terms', 'formula'  language No_ranas ~ Profundidad_charco
##   .. ..- attr(*, "variables")= language list(No_ranas, Profundidad_charco)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "No_ranas" "Profundidad_charco"
##   .. .. .. ..$ : chr "Profundidad_charco"
##   .. ..- attr(*, "term.labels")= chr "Profundidad_charco"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(No_ranas, Profundidad_charco)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "No_ranas" "Profundidad_charco"
##  $ model        :'data.frame':   20 obs. of  2 variables:
##   ..$ No_ranas          : num [1:20] 3 16 6 11 5 16 16 13 16 19 ...
##   ..$ Profundidad_charco: num [1:20] 15 39 24 19 29 49 59 36 64 49 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language No_ranas ~ Profundidad_charco
##   .. .. ..- attr(*, "variables")= language list(No_ranas, Profundidad_charco)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "No_ranas" "Profundidad_charco"
##   .. .. .. .. ..$ : chr "Profundidad_charco"
##   .. .. ..- attr(*, "term.labels")= chr "Profundidad_charco"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(No_ranas, Profundidad_charco)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "No_ranas" "Profundidad_charco"
##  - attr(*, "class")= chr "lm"

# intervalos de confianza de los parámetros
confint(modelo) 

##                         2.5 %    97.5 %
## (Intercept)        -2.2984941 5.9638854
## Profundidad_charco  0.1810911 0.3803002

rm(list = ls()) # elimina datos

set.seed(34); n <- 50; sigma <- 2; b0 <- 2; b1 <- 0.3 
x <- runif(n, 10, 30)
yhat <- b0 + b1*x
y <- rnorm(n, yhat, sd = sigma)

mod <-  lm(y~x)
par(mfrow = c(1,1))     
plot(x, y, pch = 16, las = 1, cex = 2, bty = "l", col = "skyblue", xlab = "Profundidad agua (cm)", ylab = "Número de ranas")
abline(mod, lwd = 2, col = "red")

library(arm)
nsim <- 100
bsim <- sim(mod, n.sim = nsim)
apply(coef(bsim), 2, quantile, prob=c(0.025, 0.975)) 

##       (Intercept)         x
## 2.5%    0.1600187 0.1980701
## 97.5%   3.9253881 0.3861469

quantile(bsim@sigma, prob=c(0.025, 0.975))

##     2.5%    97.5% 
## 1.701432 2.383393

quantile(coef(bsim)[,2], probs=c(0.025, 0.975))

##      2.5%     97.5% 
## 0.1980701 0.3861469

sum(coef(bsim)[,2]>1)/nsim

## [1] 0

sum(coef(bsim)[,2]>0.5)/nsim

## [1] 0

par(mfrow = c(1,2), mar=c(5,4,3,3))     
plot(x, y, pch = 16, las = 1, col = "skyblue", bty = "l", xlab = "Profundidad agua (cm)", ylab = "Número de ranas")
for(i in 1:nsim) abline(coef(bsim)[i,1], coef(bsim)[i,2], col = rgb(1, 0, 0, 0.05))

newdat <- data.frame(x = seq(10, 30, by = 0.1))
newmodmat <- model.matrix(~x, data = newdat)
fitmat <- matrix(ncol = nsim, nrow=nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- newmodmat%*%coef(bsim)[i,]

plot(x, y, pch = 16, las = 1, , col = "skyblue", bty = "l", xlab = "Profundidad agua (cm)", ylab = "Número de ranas")
abline(mod, lwd=2, col = "red")
lines(newdat$x, apply(fitmat, 1, quantile, prob=0.025), lty=3, col = "red")
lines(newdat$x, apply(fitmat, 1, quantile, prob=0.975), lty=3, col = "red")

set.seed(34) 
nsim <- 50000
bsim <- sim(mod, n.sim=nsim)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- newmodmat%*%coef(bsim)[i,]

par(mfrow = c(1,1))     
plot(x, y, pch = 16, las = 1, col = "skyblue", bty = "l", xlab = "Profundidad agua (cm)", ylab = "Número de ranas")
abline(mod, lwd=2, col = "red")
lines(newdat$x, apply(fitmat, 1, quantile, prob=0.025), lty=3, col = "red")
lines(newdat$x, apply(fitmat, 1, quantile, prob=0.975), lty=3, col = "red")

newy <- matrix(ncol=nsim, nrow=nrow(newdat)) 
for(i in 1:nsim) newy[,i] <- rnorm(nrow(newdat), mean=fitmat[,i], sd=bsim@sigma[i])
lines(newdat$x, apply(newy, 1, quantile, prob=0.025), lty=2, col = "red")
lines(newdat$x, apply(newy, 1, quantile, prob=0.975), lty=2, col = "red")

sessionInfo()

## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] arm_1.10-1    lme4_1.1-18-1 Matrix_1.2-14 MASS_7.3-50  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.0      lattice_0.20-35 digest_0.6.18   grid_3.5.1     
##  [5] nlme_3.1-137    magrittr_1.5    coda_0.19-1     evaluate_0.12  
##  [9] stringi_1.2.4   minqa_1.2.4     nloptr_1.2.1    rmarkdown_1.11 
## [13] splines_3.5.1   tools_3.5.1     stringr_1.3.1   prettydoc_0.2.1
## [17] abind_1.4-5     yaml_2.2.0      compiler_3.5.1  htmltools_0.3.6
## [21] knitr_1.20