Quiero ajustar y~intercepto+x.1*x+x.2*x^2.
Aquí está el indata. Los elementos "ab" y "ad" tienen la misma fórmula subyacente que "aa" y "ac".
library(data.table)
library(ggplot2)
set.seed(123)
DT.1 <- data.table(item= "aa", x=1:100)
DT.1[, y:=0+x*40-x^2*0.4+rnorm(100, mean=0, sd=100)]
DT.2 <- data.table(item= "ab", x=1:10)
DT.2[, y:=0+x*40-x^2*0.4+rnorm(10, mean=0, sd=100)]
DT.3 <- data.table(item= "ac", x=1:100)
DT.3[, y:=0+x*50-x^2*0.4+rnorm(100, mean=0, sd=100)]
DT.4 <- data.table(item= "ad", x=61:70)
DT.4[, y:=0+x*40-x^2*0.4+rnorm(10, mean=0, sd=100)]
DT <- rbind(DT.1, DT.2, DT.3, DT.4)
p <- ggplot(DT, aes(x=x, y=y, group=item))
p <- p + geom_point()
p <- p + facet_wrap(~item)
p
Pero al hacer un lm normal, el ajuste de "ab" y "ad" no es lo que necesito.
DT[, c("intercept", "x.1", "x.2"):=as.list(coef(lm("y~x+I(x^2)", .SD))), by = list(item)]
DT[, x.index:=1:.N, by=item]
DT.coef <- DT[x.index==1, list(item, intercept, x.1, x.2)]
DT.1 <- data.table(item="aa", x=1:100)
DT.1[, y.lm:=DT.coef[item=="aa", intercept]+x*DT.coef[item=="aa", x.1]+x^2*DT.coef[item=="aa", x.2]]
DT.2 <- data.table(item="ab", x=1:100)
DT.2[, y.lm:=DT.coef[item=="ab", intercept]+x*DT.coef[item=="ab", x.1]+x^2*DT.coef[item=="ab", x.2]]
DT.3 <- data.table(item="ac", x=1:100)
DT.3[, y.lm:=DT.coef[item=="ac", intercept]+x*DT.coef[item=="ac", x.1]+x^2*DT.coef[item=="ac", x.2]]
DT.4 <- data.table(item="ad", x=1:100)
DT.4[, y.lm:=DT.coef[item=="ad", intercept]+x*DT.coef[item=="ad", x.1]+x^2*DT.coef[item=="ad", x.2]]
DT.lm <- rbind(DT.1, DT.2, DT.3, DT.4)
p <- ggplot(DT.lm, aes(x=x, y=y.lm, group=item))
p <- p + geom_point()
p <- p + facet_wrap(~item)
p
Así que me gustaría establecer un prior que haga que los "ab" y los "ad" sean mucho más parecidos a los "aa" y los "ac".
Intentando replicar ejemplos de la documentación y viñetas de rstanarm, me quedo atascado en cómo definir el prior. Esto es lo más cerca que llego, y produce el error "Error en prior ". $location : $ no es válido para vectores atómicos".
library(rstanarm)
y.posterior <- stan_lm(y~x+I(x^2), data=DT[item=="aa", ], prior=mean(0.5), chains=1, cores=1, iter=1000, seed=12345)¿Cuál es (conceptualmente) la prioridad de y~intercepto+x.1*x+x.2*x^2? ¿Debo definir la media y la sd para el intercepto, x.1, y x.2? Puedo "dibujar" cómo debería ser la curva a priori (pico alrededor de x=50 con y alrededor de 1100), lo que se traduce en intercepto=18, x.1=40, y x.2=-0.4.
¿Cómo lo traduzco al rstanarm anterior quiere?
Espero que mi intento de utilizar métodos estadísticos avanzados con tan pocos conocimientos de ciencia estadística no insulte a nadie aquí. Por favor, no hagan comentarios sobre los datos de entrada, sólo están ahí para que el código funcione.
