Práctica 4 (2016)

Pràctica Español
Universidad Universidad Autónoma de Barcelona (UAB)
Grado Estadística Aplicada - 3º curso
Asignatura Simulació, Remostreig i Aplicacions
Año del apunte 2016
Páginas 5
Fecha de subida 27/04/2016
Descargas 3
Subido por

Vista previa del texto

Pra´ctica 4: Pruebas Permutacionales Anna Olmo, 1363013 March 13, 2016 1 Problema 1 > datos<-data.frame("Temp"=c(75.00,83.00,85.00,85.00,92.00,97.00,99.00), + "WC"=c(16.00,20.00,25.00,27.00,32.00,48.00,48.00), + "TMG"=c(1.85,1.25,1.50,1.75,1.15,1.75,1.60)) 1.1 (a) Ajusta un modelo explicando el consumo de agua (W C) a partir de la temperatura (T ) y el tiempo que se tarda en segar la hierba (T M G) por medio de un modelo lineal de regresi´on m´ ultiple.
> mod<-lm(WC~Temp+TMG,data=datos) > summary(mod) Call: lm(formula = WC ~ Temp + TMG, data = datos) Residuals: 1 2 3 4 1.0441 0.4642 -0.6935 -1.8264 5 0.1061 6 7 1.0252 -0.1197 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -121.65500 6.54035 -18.601 4.92e-05 *** Temp 1.51236 0.06077 24.886 1.55e-05 *** TMG 12.53168 1.93302 6.483 0.00292 ** --Signif. codes: 0 ✬***✬ 0.001 ✬**✬ 0.01 ✬*✬ 0.05 ✬.✬ 0.1 ✬ ✬ 1 Residual standard error: 1.245 on 4 degrees of freedom Multiple R-squared: 0.9937, Adjusted R-squared: F-statistic: 313.2 on 2 and 4 DF, p-value: 4.027e-05 1.2 0.9905 (b) Prueba la significaci´ on estad´ıstica de las variables T y T M G usando una prueba exacta de permutaciones.
> require(combinat) > Temp.mod<-summary(mod)$coefficients[2,3] > TMG.mod<-summary(mod)$coefficients[3,3] Permutamos la variable respuesta para dejar fijas las explicativas, de este modo, el consumo de agua es independiente de cuales sean sus explicativas correspondientes.
> > > > + permWC<-permn(datos$WC) newTemp<-0 newTMG<-0 for(i in 1:length(permWC)){ newmod<-lm(permWC[[i]]~Temp+TMG,data=datos) 1 + newTemp[i]<-summary(newmod)$coefficients[2,3] + newTMG[i]<-summary(newmod)$coefficients[3,3] + } > pvalorTemp<-sum(newTemp>=Temp.mod | newTemp<=-Temp.mod)/length(permWC); pvalorTemp [1] 0.001190476 > pvalorTMG<-sum(newTMG>=TMG.mod | newTMG<=-TMG.mod)/length(permWC); pvalorTMG [1] 0.009126984 > Los pvalores demuestran que ambas variables son estad´ısticamente significativas.
Otro modo de hacerlo sin usar bucles, seria mediante: > TstatLMmod<-function(x){ + mod<-summary(lm(x~Temp+TMG,data=datos)) + stat<-abs(c(mod$coefficients[2,3], + mod$coefficients[3,3])) + c(stat[1]>=Temp.mod,stat[2]>=TMG.mod) + } > apply(sapply(permWC,TstatLMmod),1,sum)/length(permWC) [1] 0.001190476 0.009126984 El resultado es el mismo que el anterior.
1.3 (c) ¿Cu´ al es la predicci´ on del consumo de agua (W C) para T = 95 y T M G = 2? Encuentra un intervalo de confianza al 95% usando el procedimiento permutacional.
Con las funciones de R seria: > # predicci´ on > predict(mod,newdata=data.frame("Temp"=95,"TMG"=2), + interval=c("prediction"),level=0.95) fit lwr upr 1 47.08295 42.41745 51.74844 > # intervalo de confianza > predict(mod,newdata=data.frame("Temp"=95,"TMG"=2), + interval=c("confidence"),level=0.95) fit lwr upr 1 47.08295 43.94917 50.21673 2 Con el proceso permutacional: > > > > + + + + + + + + > fit<-0 lower<-0 upper<-0 for(i in 1:length(permWC)){ newmod<-lm(permWC[[i]]~Temp+TMG,data=datos) fit[i]<-predict(newmod,newdata=data.frame("Temp"=95,"TMG"=2), interval=c("confidence"),level=0.95)[1] lower[i]<-predict(newmod,newdata=data.frame("Temp"=95,"TMG"=2), interval=c("confidence"),level=0.95)[2] upper[i]<-predict(newmod,newdata=data.frame("Temp"=95,"TMG"=2), interval=c("confidence"),level=0.95)[3] } cbind("fit"=mean(fit),"lwr"=mean(lower),"upr"=mean(upper)) fit lwr upr [1,] 30.85714 -0.5244408 62.23873 Otro modo de hacerlo sin usar bucles, seria mediante: > predLMmod<-function(x){ + newmod<-lm(x~Temp+TMG,data=datos) + predict(newmod,newdata=data.frame("Temp"=95,"TMG"=2), + interval=c("confidence"),level=0.95) + } > apply(sapply(permWC,predLMmod),1,mean) [1] 30.8571429 -0.5244408 62.2387265 Conclusi´ on: El m´etodo permutacional no funciona, ya que el cero est´a incluido y no deber´ıa estarlo.
2 Problema 2 > datos<-data.frame("semanas"=c(3,3,3,3,4,4,4,6,6,8,8,8,8,11,11,12,12,5,5,6,6,7,9,10,11,11,16,17,19,20,2 + "censura"=c(0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,1,1,0,1,0,1,1,1,0,1,1), + "grupo"=as.factor(c(rep("C",17),rep("N",16)))) 2.1 (a) Estudiar las dos funciones de supervivencia de manera cl´asica mediante los estimadores de Kaplan-Meier.
Para resolverlo hacer uso del paquete survival (funciones Surv, survfit y survdiff) de R. ¿Qu´e conclusiones se extraen? > > > > > require(survival) censuraR<-Surv(time=datos$semanas,event=(datos$censura-1)*(-1)) surv.mod<-survfit(censuraR~datos$grupo,type="kaplan-meier") plot(surv.mod,col=c("dark blue","dark orange"),main="",xlab="semanas",ylab="",lwd=2) legend("topright",legend=c("control","tratamiento"),fill=c("dark blue","dark orange")) 3 2.2 (b) Haz una prueba permutacional para comparar las dos fucniones de supervivencia. Esta prueba puede basarse en el test de log-rank. ¿A qu´e conclusi´ on se llega? > logRankTest<-function(dato,grupo) { + censuraR<-Surv(time=dato$semanas,event=(dato$censura-1)*(-1)) + survdiff(censuraR~grupo,rho=0)$chisq + } > LRTRUE<-logRankTest(datos,datos$grupo) Permutamos la variable grupo para dejar fijas los tiempos, de este modo, el tiempo es independiente de cual sea su grupo correspondiente > > > + + > grupo2<-t(replicate(10000,sample(datos$grupo))) newLR<-0 for(i in 1:10000){ newLR[i]<-logRankTest(datos,as.vector(grupo2[i,])) } pvalorLR<-sum(newLR>=LRTRUE)/10000; pvalorLR [1] 0.002 El pvalor demuestra que la variable es estad´ısticamente significativa.
4 ...

Tags: