Gracias Emilio por el codigo en R. Hace poco estuve revisando una situacion similar a la que ilustras en tu mensaje (i.e., seleccion de modelos e inferencia). Dale una mirada a estos dos articulos:
http://www-stat.wharton.upenn.edu/~berkr/PoSI-submit.pdf http://statweb.stanford.edu/~ckirby/brad/papers/2013ModelSelection.pdf Saludos cordiales, Jorge.- 2015-01-24 6:24 GMT+11:00 Emilio Torres Manzanera <tor...@uniovi.es>: > �Mil gracias! > Parece que el tama�o �s� que importa! aunque solo sean 4 coeficientes. > Para estimar los 4 coeficientes con cierta precisi�n, se necesitan unos > 2000 datos. La pr�xima vez que vea un estudio sobre si trabaja (s�/no) en > funci�n del sexo (hombre/mujer) y estado civil (casado/soltero), con un > tama�o de muestra de n=400, no s� qu� le dir� al autor... > > Respecto a los coeficientes, la verdad es que eran feos (eso me pasa por > copiar y pegar de otros ejemplos). Pongo abajo el c�digo definitivo que he > utilizado para que mi pobre ordenador pudiera correr la simulaci�n en un > tiempo razonable. > > �Gracias de nuevo! > Emilio > > ## ============================================================ > library(dplyr) > ## funcion auxiliar > hacertablafrecuencias <- function(x, vars=NULL, freq=NULL, sort=FALSE){ > evaldp <- function (.data, .fun.name, ..., .envir = parent.frame()) > { > args <- list(...) > args <- partial_eval(args, env = .envir) > args <- unlist(args) > if (is.function(.fun.name)) > .fun.name <- as.character(substitute(.fun.name)) > code <- paste0(.fun.name, "(.data,", paste0(args, collapse = ","), > ")") > eval(parse(text = code, srcfile = NULL)) > } > if(is.null(vars)) vars <- names(x) > nms <- c(vars, "freq") > nmsfreq <- make.names(nms, unique = TRUE)[length(nms)] > if (is.null(freq)) { > action <- paste( nmsfreq,"=n()" ) ##quote(n()) > } else { > if(freq=="freq") nmsfreq <- "freq" > vars <- vars[ vars!= "freq"] > action <- paste(nmsfreq, "=sum(",freq,")") > } > if(nmsfreq!="freq") warning(paste("freq column:",nmsfreq)) > out <- group_by_(x, .dots = vars) %>% > evaldp(summarise, action) %>% > ungroup() > if (!sort) { > out > } else { > arrange(out, desc(freq)) > } > } > > ## modelo : -1 - 4 * dat$x1 + 7 * dat$x2 - 1 *dat$x1* dat$x2 > > logisticsimulation <- function(n){ > dat <- data.frame(x1=rep(c(0:1), each=n), > x2=rep(c(0:1), time=n)) > odds <- exp(-1 - 4 * dat$x1 + 7*dat$x2 - 1 *dat$x1* dat$x2 ) > pr <- odds/(1+odds) > res <- replicate(100, { > dat$y <- rbinom(2*n,1,pr) > ttdat <- hacertablafrecuencias(dat) > coef(glm(y ~ x1*x2, data = ttdat, family = > binomial(),weights=ttdat$freq)) > }) > t(res) > } > > res <- logisticsimulation(400) # son en realidad 800 datos. Ajusta de pena > apply(res,2,median) > ## (Intercept) x1 x2 x1:x2 > ## -0.9946226 -4.2473363 19.8453136 -0.5429929 > > res <- logisticsimulation(1000) # son en realidad 2000 datos. Buen ajuste > apply(res,2,median) > ## (Intercept) x1 x2 x1:x2 > ## -1.004794 -4.120370 7.243097 -1.267561 > > > > > > > On Friday, January 23 2015, 12:23:25, Olivier Nu�ez <onu...@unex.es> > wrote: > > > Efectivamente, la normalidad tarda en llegar (un problema que merece ser > investigado) > > En cualquier caso, parece que ambos dise�os dan varianzas asint�ticas > similares de los estimadores: > > > >> n=1000 > >> > >> dat <- data.frame(x1=sample(0:1, n,replace=TRUE),x2=sample(0:1, > n,replace=TRUE)) > >> dat$intercept=1 > >> odds <- exp(-1 - 4 * dat$x1 + 7*dat$x2 - 1 *dat$x1* dat$x2 ) > >> pr <- odds/(1+odds) > >> D=diag(pr*(1-pr)) > >> X=as.matrix(dat) > >> I=t(X)%*%D%*%X > >> solve(I) > > x1 x2 intercept > > x1 0.4716741 -0.451997095 -0.014952896 > > x2 -0.4519971 0.470731667 -0.005214237 > > intercept -0.0149529 -0.005214237 0.020017370 > >> > >> dat$x1= rep(c(0,1), each = n/2) > >> dat$x2 = rep(c(0,1), times = n/2) > >> odds <- exp(-1 - 4 * dat$x1 + 7*dat$x2 - 1 * dat$x1 * dat$x2 ) > >> pr <- odds/(1+odds) > >> > >> D=diag(pr*(1-pr)) > >> X=as.matrix(dat) > >> I=t(X)%*%D%*%X > >> solve(I) > > x1 x2 intercept > > x1 0.45113285 -0.430788206 -0.014755274 > > x2 -0.43078821 0.451132851 -0.005589371 > > intercept -0.01475527 -0.005589371 0.020161833 > >> > > > > > > ----- Mensaje original ----- > > De: "Carlos J. Gil Bellosta" <c...@datanalytics.com> > > Para: "Olivier Nu�ez" <onu...@unex.es> > > CC: "Emilio Torres Manzanera" <tor...@uniovi.es>, "r-help-es" < > r-help-es@r-project.org> > > Enviados: Viernes, 23 de Enero 2015 11:14:44 > > Asunto: Re: [R-es] Simulaci�n de modelo logit con interacci�n > > > > Hola, �qu� tal? > > > > Cierto, cierto, hab�a un error en el c�digo que publiqu�. Pero el > > diagn�stico es parecido. Cuando los datos se generan con el > > coeficiente de x2 igual a 7, los coeficientes estimados tienen una > > distribuci�n extra�a, bimodal (aparentemente), en lugar de > > _normalmente_ distribuida alrededor del 7 como se espera. Supongo que > > depende del n�mero de casos en que x2 = 1 e y = 0. > > > > Un saludo, > > > > Carlos J. Gil Bellosta > > http://www.datanalytics.com > > > > El d�a 23 de enero de 2015, 9:54, Olivier Nu�ez <onu...@unex.es> > escribi�: > >> Emilio, > >> > >> espero no haberte generado mucha confusi�n con mi anterior respuesta. > >> El problema no es de separaci�n sino m�s bien de tama�o muestral. > >> Al coger el c�digo de Carlos, obten�a que "y" y "x1" eran > sistem�ticamente independiente (la tabla table(dat$y,dat$x1) tiene columnas > proporcionales). > >> En el dise�o de Carlos, el c�digo correcto para simular los datos ha de > ser dat$y= rbinom(2*n,1,pr). > >> Con tu dise�o, si elijo un tama�o muestral suficientemente grande, > obtengo estimaciones razonables de los par�metros de tu modelo: > >> > >>> res <- logisticsimulation(10000) > >>> apply(res,2,median) > >> (Intercept) x1 x2 x1:x2 > >> -0.9955877 -3.9938632 6.9967218 -0.9913839 > >> > >> Un saludo. Olivier > >> > >> > >> ----- Mensaje original ----- > >> De: "Emilio Torres Manzanera" <tor...@uniovi.es> > >> Para: r-help-es@r-project.org > >> Enviados: Jueves, 22 de Enero 2015 12:28:14 > >> Asunto: [R-es] Simulaci�n de modelo logit con interacci�n > >> > >> Hola, > >> Deseo simular un modelo logit con interacci�n, estimar sus coeficientes > y comprobar si son o no parecidos al modelo te�rico. Con este ejemplo > obtengo que los coeficientes estimados no se asemejan mucho a los > originales. �Se le ocurre a alguien cu�l es el motivo de esta discrepancia? > �y c�mo solucionarlo? > >> Muchas gracias > >> Emilio > >> > >> logisticsimulation <- function(n){ > >> dat <- data.frame(x1=sample(0:1, n,replace=TRUE), > >> x2=sample(0:1, n,replace=TRUE)) > >> odds <- exp(-1 - 4 * dat$x1 + 7*dat$x2 - 1 *dat$x1* dat$x2 ) > >> pr <- odds/(1+odds) > >> res <- replicate(100, { > >> dat$y <- rbinom(n,1,pr) > >> coef(glm(y ~ x1*x2, data = dat, family = binomial())) > >> }) > >> t(res) > >> } > >> > >> res <- logisticsimulation(100) > >> apply(res,2,median) > >> ## (Intercept) x1 x2 x1:x2 > >> ## -1.0986123 -18.4674562 20.4823593 -0.0512933 > >> > >> Deber�an salir -1, -4, 7, 1 > >> > >> _______________________________________________ > >> R-help-es mailing list > >> R-help-es@r-project.org > >> https://stat.ethz.ch/mailman/listinfo/r-help-es > >> > >> _______________________________________________ > >> R-help-es mailing list > >> R-help-es@r-project.org > >> https://stat.ethz.ch/mailman/listinfo/r-help-es > > > -- > ================================================= > Emilio Torres Manzanera > Fac. de Comercio - Universidad de Oviedo > c/ Luis Moya 261, E-33203 Gij�n (Spain) > Tel. 985 182 197 email: tor...@uniovi.es > > _______________________________________________ > R-help-es mailing list > R-help-es@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-help-es > [[alternative HTML version deleted]]
_______________________________________________ R-help-es mailing list R-help-es@r-project.org https://stat.ethz.ch/mailman/listinfo/r-help-es
