¿Cómo implementar una función de la varianza de R?

Question

Quiero estimación multivariante de la varianza de la función en R. Que es, quiero permitir la varianza (así como la media) varían de acuerdo con algún conjunto de variables independientes.

En este caso en particular, yo quiero estimar los efectos de un conjunto típico de covariables demográficas (edad, raza, educación) en la varianza de la sesión de los salarios.

Lo que es una buena manera de implementar esto en R? Hay un paquete que simplifica esto?

Puede ser que esta es sólo una búsqueda de distancia - pero después de haber buscado en el R páginas de ayuda de Google, Rseek, y StackOverflow, no puedo encontrar nada relevante en virtud de "la varianza de la función" o similar.

Cualquier sugerencias recibidas con gratitud.

---

Gracias por sus respuestas -- voy a tratar de aclarar mi pregunta.

Estoy trabajando en un máximo de probabilidad de marco. Puedo código esta de la mano de la log-verosimilitud, pero el conjunto de datos reales tiene un montón de variables y "optim" es muy lento, así que me gustaría encontrar un paquete de R que hace que esto sea más eficiente computacionalmente.

Voy a empezar con la log-verosimilitud para un básico de regresión OLS: $$ \text{ln }L = \sum (-\frac{1}{2} (\text{ln }\sigma^2 - \frac{(y - xB)^2}{\sigma^2})) $$ Entonces me relajan el supuesto de varianza constante (homoskedasticity) y redefinir la varianza como: $$ \sigma^2 = exp(Z*\gamma) $$ donde $Z$ es la matriz de las variables que afectan $\sigma^2$. (Exponentiate de modo que usted no termina con $\sigma^2$ menos que cero.) Cuando me sustituir el reajuste de parámetros de $\sigma^2$ en el original de la log-verosimilitud y el código de la nueva función de verosimilitud logarítmica en R, me sale esto:

ll.normal.vary <- function (par, X, Y, Z) {
  beta  <-par[1:ncol(X)]
  gamma <- par[(ncol(X)+1):(ncol(X)+ncol(Z))]   
  -1/2* sum((Z %*% gamma) + ((Y - X %*% beta)^2)/exp(Z %*% gamma))
}

Entonces me optimizar:

v.optim1 <- optim (par = start1, fn=ll.normal.vary, X=x.mat, Y=y.vec, Z=z.mat, 
                   method = "BFGS", hessian = F, control = list(fnscale = -1))
v.optim1$par
v.optim1$value

Aquí están algunos datos de ejemplo si quieres probar:

var1   <- c(0,0,0,1,1,0)
var2   <- c(.28, .07, -.05, .38, .08, -.1)
var3   <- c(-.11, -.17, -.17, -.05, .1, -.01)
x.mat  <- cbind(var1, var2, var3)
y.vec  <- c(.46, .77, .49, .59, .60, .44)
z.mat  <- cbind(var1, var2) 
start1 <- rep(0.1, ncol(x.mat)+ncol(z.mat))

Gracias de nuevo por los consejos.

Answer 1

De forma iterativa re-mínimos cuadrados ponderados es probablemente la respuesta. Comience con un conjunto de pesas todos iguales a 1, o igual a la inversa de su mejor estimación de la varianza en cada punto. El uso de mínimos cuadrados ponderados para ajustar el modelo con el que el conjunto de pesos. A continuación, vuelva a calcular su peso basado en el último set de parámetros; y reformar su modelo con ese conjunto de pesos. Continuar hasta lograr la convergencia. Esto es, básicamente, el método comúnmente utilizado para el ajuste de modelos lineales generalizados, donde la varianza en cualquier punto depende de la función de enlace y la distribución de la respuesta.

Sin duda hay funciones R que hacer esto para usted, pero si usted quiere jugar con diferentes modelos de varianza se relaciona con sus diferentes parámetros que usted puede encontrar que es igual de fácil hacerlo usted mismo.

Alternativamente, si usted tiene un bastante estándar, la varianza del modelo, usted puede encontrar que sólo el glm() función, con una adecuada función de enlace y de distribución, que funciona para usted. El "cuasi" la familia puede ser particularmente útil, ya que permite especificar tanto la función de enlace y la relación de la varianza de la respuesta, con un mínimo de supuestos acerca de los otros momentos de la respuesta.