Simulacro_Porras

Librerias necesarias

library(tidyverse)

Lectura de datos

casoEF <- readxl::read_excel(path = 'Caso EF.xlsx')

datos <- 
  casoEF %>% 
  dplyr::mutate(genero = dplyr::recode(genero, 'femenino' = 0, 'masculino' = 1))

Pregunta 01

Use el conjunto de datos casoEF.csv.

Se desea utilizar un Analisis de Regresion Lineal Multiple Utilice solo los datos de las personas pertenecientes a la región A y la(s) variable(s) predictora(s) significativa(s) a un nivel de significacion de 0.01.

Antes de correr su funcion, donde crea conveniente fije: RNGkind(sample.kind=“Rounding”) set.seed(55) En todos los casos presente sus resultados con aproximacion a 3 decimales.

Filtro de datos, region A

datosA <- 
  datos %>% 
  dplyr::filter( region %in% c('A') ) %>% 
  dplyr::select(-c(region))

Valor de semilla

semilla <- 55

El valor del error de prediccion es

lm01 <- lm(formula = peso ~ ., data = datosA)
summary(lm01)

## 
## Call:
## lm(formula = peso ~ ., data = datosA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6925 -0.1955  0.0163  0.1636  0.8128 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.242e+02  2.181e+01  -5.696 6.89e-08 ***
## genero      -2.001e-02  5.022e-02  -0.398    0.691    
## edad         7.895e-04  3.072e-03   0.257    0.798    
## calorias     8.558e-02  1.504e-02   5.689 7.11e-08 ***
## estatura     9.369e-02  1.015e-01   0.923    0.358    
## ejercicios   1.600e-01  3.610e-02   4.432 1.86e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3019 on 141 degrees of freedom
## Multiple R-squared:  0.9989, Adjusted R-squared:  0.9989 
## F-statistic: 2.54e+04 on 5 and 141 DF,  p-value: < 2.2e-16

data.sig <- 
  datosA %>% 
  dplyr::select(c(peso, calorias, ejercicios))

# Variables significativas: Calorias y ejercicios
lm02 <- lm(formula = peso ~ ., data = data.sig)
VEP <- mean( (lm02$residuals)^2 );VEP

## [1] 0.088063

Error de prediccion usando PRESS

PRESSp<-function(datos,y){
  datos<-as.matrix(datos)
  n<-nrow(datos)
  resi<-lm(datos[,y]~datos[,-y])$res
  APE<-sum(resi^2)/n
  resid<-c()
  for(i in 1:n){
    betase<-lm(datos[-i,y]~datos[-i,-y])$coe
    xi<-c(1,datos[i,-y])
    yest<-sum(betase*xi)
    resid[i]<-datos[i,y]-yest
  }
  pressp<-sum(resid^2)/n
  return(list(APE=APE,pressp=pressp))
}

RNGkind(sample.kind='Rounding')
set.seed(semilla)
PRESSp(datos = data.sig, y = 1)

## $APE
## [1] 0.088063
## 
## $pressp
## [1] 0.09103554

Error de prediccion usando validacion cruzada

crossval<-function(datos,K,r,d){
  datos<-as.matrix(datos)
  n<-nrow(datos)
  EVC<-c()
  resid<-c()
  subm<-floor(n/K)
  resi<-lm(datos[,d]~datos[,-d])$res
  APE<-sum(resi^2)/n
  for(i in 1:r){
    indices<-sample(n,n)
    azar<-datos[indices,]
    
    for(j in 1:K){
      unid<- ((j-1)*subm+1):(subm*j)
      if (j==K)
      {
        unid<-((j-1)*subm+1):n
      }
     datosp<-azar[unid,]
     datose<-azar[-unid,]
     ye<-datose[,d]
     xe<-datose[,-d]
     betas<-lm(ye~xe)$coef
     datosp1<-cbind(1,datosp[,-d])
     estim<-datosp1%*%betas
     resid[j]<-sum((datosp[,d]-estim)^2)
    }
    EVC[i]<-sum(resid)/n
  }
  EVCP<-mean(EVC)
  cvEVC<-sd(EVC)*100/EVCP
  sesgo<-EVCP-APE
  return(list(APE=APE,EVCP=EVCP,EVCsd = sd(EVC), cvEVC=cvEVC,sesgo=sesgo))
}

RNGkind(sample.kind="Rounding")
set.seed(semilla)
crossval(datos = data.sig, K = 7, r = 27, d = 1)

## $APE
## [1] 0.088063
## 
## $EVCP
## [1] 0.09130565
## 
## $EVCsd
## [1] 0.0009864908
## 
## $cvEVC
## [1] 1.080427
## 
## $sesgo
## [1] 0.003242647

kopt<-function(datos,K,r,d){
  resul<-c()
  t<-length(K)
  for(i in 1:t){
    resul[i]<-system.time(crossval(datos,K[i],r,d))[3]
  }
  plot(K,resul,type="b",col=4)
  return(resul)
}

ks<-c(3,5,7,10,12,14)
RNGkind(sample.kind="Rounding")
set.seed(semilla)
kopt(data.sig,ks,80,1)

## [1] 0.11 0.19 0.26 0.37 0.43 0.51

Error de prediccion con validacion cruzada generalizada

vcg<-function(datos,y){
  datos<-as.matrix(datos)
  n<-nrow(datos)
  modelo<-lm(datos[,y]~datos[,-y])
  resi<-modelo$res
  APE<-sum(resi^2)/n
  Hp<-hatvalues(modelo)
  vcg<-APE/(1-sum(Hp)/n)^2
  return(list(APE=APE,vcg=vcg))
}

vcg(datos = data.sig, y = 1)

## $APE
## [1] 0.088063
## 
## $vcg
## [1] 0.09177051

Error de prediccion con Bootstrap simple

bs<-function(datos,B,y){
  datos<-as.matrix(datos)
  n<-nrow(datos)
  modelo<-lm(datos[,y]~datos[,-y])
  resi<-modelo$res
  APE<-sum(resi^2)/n
  bootape<-c()
  for(i in 1:B){
    indices<-sample(n,n,T)
    modelob<-lm(datos[indices,y]~datos[indices,-y])$coe
    yest<-cbind(1,datos[,-y])%*%modelob
    bootape[i]<-sum((yest-datos[,y])^2)/n
  }
  EPBS<-mean(bootape)
  return(list(APE=APE,EPBS=EPBS))
}

RNGkind(sample.kind="Rounding")
set.seed(semilla)
bs(datos = data.sig, B = 148, y = 1)

## $APE
## [1] 0.088063
## 
## $EPBS
## [1] 0.08948012

Error de prediccion con Bootstrap refinado

br<-function(datos,B,y){
  datos<-as.matrix(datos)
  n<-nrow(datos)
  modelo<-lm(datos[,y]~datos[,-y])
  resi<-modelo$res
  APE<-sum(resi^2)/n
  bootape1<-c()
  bootape2<-c()
  for(i in 1:B){
    indices<-sample(n,n,T)
    modelob<-lm(datos[indices,y]~datos[indices,-y])
    coeboot<-modelob$coe
    resboot<-modelob$res
    yest<-cbind(1,datos[,-y])%*%coeboot
    bootape1[i]<-sum((yest-datos[,y])^2)/n
    bootape2[i]<-sum(resboot^2)/n
  }
  APE1<-mean(bootape1)
  APE2<-mean(bootape2)
  BR<-APE+APE1-APE2
  return(list(APE=APE,EPBS=APE1,EPBR=BR))
}

RNGkind(sample.kind = "Rounding")
set.seed(semilla)
br(datos = data.sig, B = 148, y = 1)

## $APE
## [1] 0.088063
## 
## $EPBS
## [1] 0.08948012
## 
## $EPBR
## [1] 0.09050549

Error de prediccion mediante Bootstrap

boot0.632<-function(data,B,d)
{
  data<-as.matrix(data)
  n<-dim(data)[1]
  p<-dim(data)[2]
  y<-data[,d]
  x<-data[,-d]
  reg<-lm(y~x)
  APE<-sum(reg$resi^2)/n
  one<-rep(1,n)
  x1<-cbind(one,x)
  betas<-matrix(0,B,p)
  indices<-matrix(0,B,n)
  resid2b<-rep(0,n)
  for(i in 1:B){
    ind<-sample(1:n,n,replace=TRUE)
    indices[i,]<-ind
    betas[i,]<-lm(y[ind]~x[ind,])$coef
  }
  for(r in 1:n){
    c1<-0
    indices2<-matrix(0,B,n)
    ID=rep(0,B)
    for(i in 1:B){
      for(j in 1:n){
        if(indices[i,j]==r) indices2[i,j]=1
      }
    }
    ind<-apply(indices2,1,sum)
    for(i in 1:B) if(ind[i]==0) ID[i]=1
    n1<-sum(ID)
    v<-rep(0,n1)
    for(i in 1:B){
      if(ID[i]==1) {
        c1=c1+1
        v[c1]=i
      }
    }
    resi2<-matrix(0,n1,n)
    for(u in 1:n1){
      betas1<-betas[v[u],]
      yp<-x1%*%betas1
      resi<-(y-yp)
      resi2[u,]<-resi^2
    }
    residb<-apply(resi2,1,sum)
    resid2b[r]<-mean(residb)
  }
  e0<-mean(resid2b)/n
  e0.632<-0.368*APE+0.632*e0
  return(list(APE=APE,e0=e0,e0.632=e0.632))
}

RNGkind(sample.kind ="Rounding")
set.seed(semilla)
boot0.632(data = data.sig, B = 148, d = 1)

## $APE
## [1] 0.088063
## 
## $e0
## [1] 0.08948808
## 
## $e0.632
## [1] 0.08896365

Pregunta 02

Considere solo los datos de la región A.

Se desea aplicar un Analisis de Regresion Lineal Simple para explicar el peso de una persona en funcion de la variable predictora edad. Considere como estadistico de interés a la pendiente de la recta.

Antes de correr su funcion, donde crea conveniente fije: RNGkind(sample.kind=“Rounding”) set.seed(89)

En todos los casos presente sus resultados con aproximacion a 3 decimales.

datos02 <- 
  dplyr::select(datosA, c(peso, edad))

Semilla

semilla02 <- 89

# Jackknife
jack2.2<-function(datos,estadistico,...){
  n<-nrow(datos)
  estjack<-c()
  for(i in 1:n){
    estjack[i]<-estadistico(datos[-i,],...)
  }
  esjack<-mean(estjack)
  eejack<-(n-1)*sd(estjack)/sqrt(n)
  resu<-c(esjack, eejack)
  return(resu) 
}

# Intervalo de confianza Jackknife
jack2.2.ic<-function(datos,B,conf,estadistico,...){
  alfa<-1-conf
  n<-nrow(datos)
  estjack1<-matrix(0,B,2)
    for(i in 1:B){
    indices<-sample(n,n,T)
    mboot<-datos[indices,]
    estjack1[i,]<-jack2.2(mboot,estadistico,...)
    }
  esta<-estadistico(datos,...)
  eejackm<-jack2.2(datos,estadistico,...)[2]
  pivot<-(estjack1[,1]-esta)/estjack1[,2]
  LI.me<-esta-qnorm(1-alfa/2)*eejackm
  LS.me<-esta+qnorm(1-alfa/2)*eejackm
  ME<-c(LI.me,LS.me)
  LI.mes<-esta+quantile(pivot,alfa/2)*eejackm
  LS.mes<-esta+quantile(pivot,1-alfa/2)*eejackm
  MES<-c(LI.mes,LS.mes)
  return(list(ME=ME,MES=MES))
}

lm.pend <- function(datos) {
  
  datos.matrix <- as.matrix(datos)
  pendiente <- lm(datos.matrix[, 1] ~ datos.matrix[, -1])$coefficients[2]
  return(pendiente)
  
}

RNGkind(sample.kind = "Rounding")
set.seed(semilla02)
jack2.2.ic(datos02, B = 148, conf = 0.94, lm.pend)

## $ME
## datos.matrix[, -1] datos.matrix[, -1] 
##         -0.1734490          0.1671998 
## 
## $MES
## datos.matrix[, -1] datos.matrix[, -1] 
##         -0.1994212          0.1621976

# system.time(jack2.2.ic(datos02, B = 148, conf = 0.94, lm.pend))[3]

BCA

bca <- function(datos,B,conf,estadistico,...){
  alfa<-1-conf
  n<-nrow(datos)
  estb<-c()
  estjack<-c()
  esta<-estadistico(datos,...)
  for(i in 1:B){
   estb[i]<-estadistico(dplyr::sample_n(tbl = datos, size = n, replace = T), ...)
  }
  z0 <- qnorm(mean(estb < esta))
  
  for(j in 1:n){
    estjack[j]<-estadistico(datos[-j, ],...)
  }
  
  num<-sum((mean(estjack)-estjack)^3)
  den<-6*((sum((mean(estjack)-estjack)^2))^1.5)
  a<-num/den
  alfa1<-pnorm(z0+(z0+qnorm(alfa/2))/(1-a*(z0+qnorm(alfa/2))))
  alfa2<-pnorm(z0+(z0+qnorm(1-alfa/2))/(1-a*(z0+qnorm(1-alfa/2))))
  LI<-quantile(estb,alfa1)
  LS<-quantile(estb,alfa2)
  lim<-c(LI,LS)
  return(lim)
}


RNGkind(sample.kind = "Rounding")
set.seed(semilla02)
bca(datos = datos02, B = 148, conf = 0.94, estadistico = lm.pend)

##  6.302328%  98.73415% 
## -0.1340222  0.1932282

Jacknife despues de bootstrap

boot2<-function(datos,B,estadistico,...){
  n<-nrow(datos)
  estboot<-c()
  for(i in 1:B){
    indices<-sample(n,n,T)
    estboot[i]<-estadistico(datos[indices,],...)
  }
  esta<-estadistico(datos,...)
  esboot<-mean(estboot)
  eeboot<-sd(estboot)
  seboot<-(esboot-esta)
  return(list(esboot=esboot,eeboot=eeboot,seboot=seboot)) 
}

jdb<-function(datos,B,estadistico,...){
  n<-nrow(datos)
  medi<-c()
  for(i in 1:n){
    medi[i]<-boot2(datos[-i,],B,estadistico,...)$seboot
  }
  vari<-(n-1)^2*var(medi)/n
  return(vari)
}

RNGkind(sample.kind = "Rounding")
set.seed(semilla02)
jdb(datos = datos02, B = 148, estadistico = lm.pend)

## [1] 0.008337628