Ejercicio Muestreo por lotes (uniparamétrico). Repaso

Se seleccionan aleatoriamente 1000 muestras de suelo que se combinan en lotes de tamagno 10 muestras. Los lotes son examinados para evaluar la presencia de una toxina, es decir, las muestras individuales de suelo no se examinan. Concretamente, en 100 lotes, de 10 muestras de suelo cada uno, se detecto la toxina en 12 lotes.

  1. Calcular el EMV para \(p=\text{probabilidad de que la toxina esté presente en una muestra de suelo individual}.\)

¿Qué datos observamos? ¿Cuál es el modelo estadístico que subyace en estos datos? \(Y_1,\dots,Y_{100}\sim Modelo(\theta)\)

¿Cómo se relaciona \(\theta\) con \(p\)? ¿Si se conoce \(p\), cuanto vale \(\theta\) (o al revés)? ¿\(P(Y_1=1)\) o \(P(Y_1=0)\)?

# EMV y de theta

sY<-12
n<-100

mlv<-function(theta){
  return(-sY*log(theta)-(100-sY)*log(1-theta))
}
minimoTheta<-optim(0.5,mlv,method="Brent",lower=0.01,upper=0.99,hessian=TRUE)
emvTheta<-minimoTheta$par
emvTheta #sY/n
## [1] 0.12
# EMV de p 
g<-function(theta){
  return(1-(1-theta)^0.1)
}

#p<-=g(theta)=1-(1-theta)^(1/10) es invariante por transformacion
emvP<-g(emvTheta)
emvP
## [1] 0.01270198
  1. Determinar la distribución asintótica del EMV de \(p\) y calcular IC aproximados de Wald para \(p\) con confianza 0.95.
# Informacion de Fisher observada
infObsTheta<-minimoTheta$hessian
infObsTheta # 1/((emvTheta*(1-emvTheta))/n)
##          [,1]
## [1,] 947.0857
# Delta metodo

gPrima<-function(theta){
  return(0.1*(1-theta)^-0.9)
}

emvVarP<-(1/infObsTheta)*gPrima(emvTheta)^2
emvVarP
##              [,1]
## [1,] 1.329052e-05
#IC Wald 95% no son invariantes por transformacion
seP<-sqrt(emvVarP)
emvP-qnorm(0.975)*seP
##             [,1]
## [1,] 0.005556701
emvP+qnorm(0.975)*seP
##            [,1]
## [1,] 0.01984725
  1. Calcular ICRV con confianza aproximada de 0.95 para \(p\).
#IC RV para theta

lv<-function(theta){
  return(sY*log(theta)+(100-sY)*log(1-theta))
}
d<--minimoTheta$value-qchisq(0.95,1)/2
d
## [1] -38.61323
# ICs para theta
plot(seq(0.001,0.4,length.out=5000),lv(seq(0.001,0.4,length.out=5000)),type="l")
abline(h=d)
abline(v=0.12,col=2) # emv
abline(v=c(0.066,0.193),lty=1) # IC RV

#IC RV para p 95% son invariantes por transformacion
g(0.066);g(0.193)
## [1] 0.006804627
## [1] 0.02121489
# Comparacion de los intervalos. Ventajas e inconvenientes
# c(0.0056,0.0198),lty=2) # IC Wald, siempre simetrico 
# c(0.0068,0.0212),lty=1) # IC RV

Ejercico genética de Fisher. Adaptación Ejercicio 9 (multiparamétrico).

Se cuenta con una muestra de tamaño n=3839 de una población que toma 4 valores (fenotipos ab, Ab, aB, AB) distintos con probabilidades \(p_1,\) \(p_2,\) \(p_3,\) y \(p_4\) tales que \(\sum_{r=1}^4 p_r=1\). Concretamente se observan las siguientes frecuencias para las tres primeras clases: 1997, 904 y 906. La hipótesis de Fisher (en contra de lo propuesto por Mendel) es que las probabilidades dependen de un solo parámetros \(\theta\) de la forma siguiente: \[H_0: p_1=\frac{2+\theta}{4}, p_2=\frac{1-\theta}{4}=p_3, p_4=\frac{\theta}{4}.\] a. Calcular el EMV para \(p=(p_1,p_2,p_3)^{'}\) e IC de Wald con confianza aproximada del 95% para \(p_1\), \(p_2\) y\(p_3.\)

n<-3839
fr<-c(1997,904,906)
mlv<-function(p){
  return(-sum(fr*log(p))-(n-sum(fr))*log(1-sum(p)))
}

# EMV para p con optim y Nelder-Mead (por defecto)

# optim minimiza -logver ~ maximiza -logver. Argumentos:
# par: valores iniciales de los parametros
# fn: funcion a minimizar
# hessian=TRUE calcula la matriz de derivadas segundas de fn ene l EMV
# method: por defecto en el caso multiparamétrico Nelder-Mead. Uniparamétrico: Brent (lower,upper)
vi<-rep(0.3,3)
minimo<-optim(par=vi,fn=mlv,hessian=TRUE)
## Warning in log(1 - sum(p)): Se han producido NaNs

## Warning in log(1 - sum(p)): Se han producido NaNs

## Warning in log(1 - sum(p)): Se han producido NaNs

## Warning in log(1 - sum(p)): Se han producido NaNs

## Warning in log(1 - sum(p)): Se han producido NaNs

## Warning in log(1 - sum(p)): Se han producido NaNs

## Warning in log(1 - sum(p)): Se han producido NaNs

## Warning in log(1 - sum(p)): Se han producido NaNs
# EMV 

emv<-minimo$par# EMV
emv
## [1] 0.5202219 0.2354914 0.2359584
p1Hat<-fr[1]/n;p2Hat<-fr[2]/n;p3Hat<-fr[3]/n;p4Hat<-1-p1Hat-p2Hat-p3Hat
p1Hat;p2Hat;p3Hat
## [1] 0.5201875
## [1] 0.235478
## [1] 0.235999
# Estimador de la varianza

infObs<-minimo$hessian # Matriz IF observada. OJO: menos Hesianno ~ mlv. 3x3
emvVar<-solve(infObs) # Inversa de IF observada. Caso uniparametrico 1/infObs
emvVar
##               [,1]          [,2]          [,3]
## [1,]  6.499573e-05 -3.192230e-05 -3.197831e-05
## [2,] -3.192230e-05  4.689324e-05 -1.447523e-05
## [3,] -3.197831e-05 -1.447523e-05  4.695013e-05
# Error estandar

se<-sqrt(diag(emvVar)) # Errores estandar de EMV 
se
## [1] 0.008061993 0.006847864 0.006852016
# IC para p1, p2 y p3
tabla.Wald<-cbind(emv,se,extremo.inferior=emv-qnorm(0.975)*se,
    extremo.superior=emv+qnorm(0.975)*se)
nombres.par<-c("p1","p2", "p3")
rownames(tabla.Wald)<-nombres.par
round(tabla.Wald,4) #usar 4 decimales
##       emv     se extremo.inferior extremo.superior
## p1 0.5202 0.0081           0.5044           0.5360
## p2 0.2355 0.0068           0.2221           0.2489
## p3 0.2360 0.0069           0.2225           0.2494
  1. Calcular p-valor del test de Wald para contrastar \(H_0\).
# Hipótesis compuesta: H0: p2-p3=0; p1+p2-3/4=0
#g1(p)=p2-p3=0
#g2(p)=p1+p2-3/4=0

g1<-function(p){
  p1<-p[1];p2<-p[2];p3<-p[3]
  return(p2-p3)
}

g2<-function(p){
  p1<-p[1];p2<-p[2];p3<-p[3]
  return(p1+p2-3/4)
}


# Delta método
G<-matrix(c(0,1,-1,1,1,0),2,3,byrow=TRUE)
varG<-G%*%emvVar%*%t(G)
Wobs<-cbind(g1(emv),g2(emv))%*%solve(varG)%*%t(cbind(g1(emv),g2(emv)))
Wobs
##          [,1]
## [1,] 2.043841
pvalW<-1-pchisq(Wobs,2)
pvalW # No se rechaza H0
##           [,1]
## [1,] 0.3599031
  1. Calcular IC para \(\theta\) con confianza aproximada de \(0.95\) basado en el estadístico de Wald.
# Bajo H0 se trata IC Wald para theta uniparamétrico: EMV de theta bajo H0 y su distribucion asintótica

#mlv bajo H0, solo depende de theta
mlvH0<-function(theta){
  pH0<-c((2+theta)/4,(1-theta)/4,(1-theta)/4)
  return(-sum(fr*log(pH0))-(n-sum(fr))*log(1-sum(pH0)))
}
viH0<-0.1
minimoH0<-optim(viH0,mlvH0,method="Brent",lower=0.01,upper=0.99,hessian=TRUE)
emvH0<-minimoH0$par
emvH0
## [1] 0.0357123
infObsH0<-minimoH0$hessian
emvVarH0<-1/infObsH0

emvH0-qnorm(0.975)*sqrt(emvVarH0)
##            [,1]
## [1,] 0.02390586
emvH0+qnorm(0.975)*sqrt(emvVarH0)
##            [,1]
## [1,] 0.04751874
  1. Calcular p-valor del test de RV para contrastar \(H_0\).
Qobs<-2*(-minimo$value+minimoH0$value)
Qobs
## [1] 2.018659
pvalQ<-1-pchisq(Qobs,2)
pvalQ # Misma conclusion pvalor Wald, aunque no tienen por que coincidir
## [1] 0.3644632
  1. Calcular IC para \(\theta\) con confianza aproximada de \(0.95\) basado en el estadístico de RV.
lvH0<-function(theta){
  pH0<-c((2+theta)/4,(1-theta)/4,(1-theta)/4)
  return(sum(fr*log(pH0))+(n-sum(fr))*log(1-sum(pH0)))
}
d<--minimoH0$value-qchisq(0.95,1)/2
plot(0,0,type="n",xlim=c(0,0.1),ylim=c(-4070,-4080))
ss<-seq(0,0.1,length.out=10000)
for(i in 1:10000){
points(ss[i],lvH0(ss[i]),pch=16,lwd=0.5)
}
abline(h=d)
abline(v=c(0.025,0.049),col=2)

Ejercicio propuesto

Sean \(X_1,\dots,X_{50}\) v.a.i.i.d. con distribución discreta que toma los valores 1, 2, 3, 4 y 5 con probabilidades \(p_1,\) \(p_2,\) \(p_3,\) \(p_4\) y \(p_5\) respectivamente, con \(p_r\geq 0, r=1,\dots,5\) y \(\sum_{r=1}^5p_r=1.\) Sea \(Y_r=\sum_{i=1}^{50}\mathbb{I}_{(X_i=r)},\) \(r=1,\dots,5.\) Si se observa \(Y_1=7, Y_2=6, Y_3=8, Y_4=16\) e \(Y_5=13.\)

  1. Obtener el p-valor del test de Wald para contrastar la hipótesis \(H_0: p_1=p_2=p_3, p_4=p_5\) contra \(H_1:\) el resto.
  2. Obtener el p-valor del test de la razón de verosimilitudes para contrastar la hipótesis \(H_0: p_1=p_2=p_3, p_4=p_5\) contra \(H_1:\) el resto.
  3. Obtener el p-valor del test de Wald para contrastar la hipótesis \(H_0: p_2-2p_1=0, p_3=p_4\) contra \(H_1:\) el resto.
  4. Obtener el p-valor del test de la razón de verosimilitudes para contrastar la hipótesis \(H_0: p_2-2p_1=0, p_3=p_4\) contra \(H_1:\) el resto.