Principios y Aplicaciones

Mate científicoMate científico "At least once in a lifetime it is convenient to put everything to discussion" (Descartes)

Readmore

useRLengua franca Usted no puede, no debe dejar de tenerlo...

Readmore

SuperCarlitosSuperCarlitos Y pa mi que la evolución tuvo que empesá má o meno así...

Readmore

Libros para descargar (gratis) sobre Diseño y Análisis Experimental

sábado, 27 de noviembre de 2010

0
Logit con variables explicativas cuantitativas. SPSS

0
Logit con variable respuesta multinomial/ordinal. SPSS

0
Logit con variables explicativas categóricas. SPSS

jueves, 25 de noviembre de 2010

0
Análisis de datos funcionales (FDA). Parte IV: regresión PLS

0
Análisis de datos funcionales (FDA). Parte III: regresión funcional

0
Análisis de datos funcionales (FDA). Parte II: Aproximación de curvas

0
Análisis de datos funcionales (FDA). Parte I: Datos funcionales

0
Diseño factorial. Aplicación en R.

0
Diseño completamente aleatorizado por bloques en R

Blog_DISEÑOS DE BLOQUES COMPLETOS e INCOMPLETOS

Ejemplo en R.
 Block=rep(c("P15","P20","P25","P30","P35"),c(4,4,4,4,4))

domingo, 13 de diciembre de 2009

0
Modelos de regresión logística o logit

Modelos logit
con variables explicativas cuantitativas observadas sin error

 


 I.            Formulación del modelo


·         Modelo de regresión logística simple
Consideraremos el caso de una única variable explicativa cuantitativa X para una variable aleatoria binaria Y. Utilizaremos un modelo lineal para el logaritmo de la ventaja de respuesta Y=1 en cada valor observado de x de la variable explicativa:

Ln[p(x)/(1-p(x))] = alfa + beta*x

Que equivalentemente se puede expresar de la siguiente forma en términos de probabilidad de respuesta 1 en x:

p(x)= e(alfa+beta*x)/(1+ e(alfa+beta*x))

curva de respuesta que es estrictamente creciente si beta>0 y estrictamente decreciente para beta <0.

Interpretación de sus parámetros
a.       Para beta=0,
          i.      tenemos p(x)= e(alfa)/(1+ e(alfa)), es decir, la variable Y es independiente de X.
          ii.      alfa es el valor común del logaritmo de las ventajas de respuesta Y=1 frente a la respuesta Y=0
b.      alfa se puede interpretar como el valor del logaritmo de la ventaja de respuesta 1 para un individuo con X=0
c.       por cada unidad de incremento en X, el logit de respuesta 1 aumenta aditivamente beta unidades. La ventaja de la respuesta 1 en cada x observado es: p(x)/(1-p(x)) = e(alfa+beta*x) = e(alfa)*(e(beta)^x), así la ventaja de la respuesta 1 aumenta multiplicativamente e(beta) por cada unidad de incremento en X.
d.      Cociente de ventajas de respuesta 1 para dos valores diferentes x1 y x2 de X: theta(x1,x2)= [ p(x1)/(1+p(x1)) ]/[ p(x2)/(1+p(x2)) ] = e(beta*(x1-x2)), que pertenece a (0, infinito) y:
          i.       theta(x1,x2)= 1 sii p(x1)=p(x2).
        ii.      theta(x1,x2)> 1 sii p(x1)>p(x2), la ventaja de respuesta 1 es e(beta*(x1-x2)) veces mayor para X=x1 que para X=x2.
        iii.      theta(x1,x2)< 1 sii p(x1)la ventaja de respuesta 1 es 1/e(beta*(x1-x2)) veces mayor para X=x2 que para X=x1.

Y para dos valores de X que se diferencian en una unidad: theta(x+1,x)= theta(deltaX=1) = e(beta)


·         Modelo de regresión logística múltiple
Si consideramos el caso de R variables explicativas cuantitaticas no aleatorias (X1, X2, …, XR), para cada combinación de valores observados X1=x1, X2=x2, …, XR=xR de las variables explicativas, la variable respuesta Y tiene distribución de Bernoulli, con p(x1,…,xR) = P[Y=1|X1=x1,…,XR=xR]=E[Y|X1=x1,…,XR=xR]. Entonces la fórmula del modelo es la siguiente:
Y(x1,…,xR)=p(x1,…,xR)+épsilon(x1,…,xR) 
donde épsilon(x1,…,xR) son errores aleatorios que se consideran centrados e independientes. Así, tenemos:
p(x1,…,xR)= e(alfa+sum(beta_r*x_r))/(1+exp(alfa+sum(beta_r*x_r)))

Si llamamos alfa=beta_0, X=(X0,X1,…,XR)` y x=(x0,x1,…,xR)`, beta=(beta0,beta1,…,betaR)’  con  X0=1, tenemos:
p(x)= e(sum(beta_r*x_r))/(1+e(sem(beta_r*x_r))) = e(beta’*x)/(1+e(beta’*x))
o equivalentemente ln[p(x)/(1-p(x))] = sum(beta_r*x_r)

Interpretación de sus parámetros

a.       Para beta_r=0, para todo r=1,…,R,
          i.      tenemos p(x)=e(beta_0)/(1+e(beta_0)), Y es independiente de las variables explicativas.
         ii.      beta_0 es el valor común del logaritmo de las ventajas de respuesta =1 frente a Y=0.
b.      Cociente de ventajas de respuesta Y=1 para dos combinaciones diferentes de valores de las variables explicativas x1=(1,x11,…,x1R)’ y x2=(1,x21,…,x2R)’ es: theta(x1,x2) = [p(x1)/(1-p(x1))]/[p(x2)/(1-p(x2))] = e(sum(beta_R(x1r-x2r)))

Y para dos valores de X que se diferencian en una unidad: theta(x+1,x)=theta(deltaX1=1,…,deltaXR=1)=e(sum(beta_r))=prod(e(beta_R)),aal aumentar una unidad de una variable y controlar las demás, la ventaja de respuesta 1 queda multiplicada por la exponencial del coeficiente de la variable incrementada.
          i.      Exponencial de un parámetro>1, entonces la probabilidad de respuesta =1 aumenta cuando aumenta la variable correspondiente y se controlan las demás
           ii.      Exponencial de un parámetro<1, se cumple la relación inversa.
   

·         Modelos con interacción
La posibilidad de interacción entre las variables explicativas de un modelo de regresión logística múltiples implica que lso cocientes de ventajas que miden la asociación entre la variable de respuesta y cada variable explicativa ya no son independientes del valor fijo del resto de variables explicativas controladas. Esto significa que los modelos anteriormente analizados son modelo sin interacción porque el grado de asociación entre la variable de respuesta y cada una de las variables explicativas es el miso en todas las combinaciones de niveles de las otras variables independientes.

Existen interacciones de distintos órdenes:
         i.            Orden uno (entre dos variables explicativas): la asociación entre la variable de respuesta y una variable, depende de los valores de una tercera que interacciona con ésta última. 
         ii.            Orden dos: involucran a tres variables

La interacción entre dos variables cuantitativas se incluye en el modelo de regresión logística múltiple como producto de ambas variables.
Ln[p(x)/(1-p(x))] = sum(beta_r*x_r) + sum(sum(beta_rs*x_r*x_s))

 Interpretación de sus parámetros
a.       El término de interacción entre dos variables cuantitativas Xr y Xs es de la forma beta_rs*Xr*Xs.
b.      El Cociente de ventajas de respuesta =1 cuando se incrementa en una unidad una variable y se controlan fijas las demás, depende del valor de las variables controladas: theta(deltaXl=1|x1,…,xl-1,xl+1,…,xR) = e(beta_l+sum(beta_lr*x_r)), para todo l=1,…,R.

 Nomenclatura:
  • Variable de confusión: está asociada con el factor de riesgo de modo que la asociación marginal entre la variable de respuesta y el factor de riesgo cambia significativamente al incluirla en el análisis estadístico. Tienen que ser considerados en el modelo aunque pueden no interaccionar con el factor de riesgo.
  • Variable modificadora: modifica el efecto cuando la asociación entre la variable de respuesta y el factor de riesgo cambia en función de sus valores. Es una variable que interacciona con el factor de riesgo.

II.            Ajuste
Disponemos de N observaciones (tamaño muestral) de N variables de Bernoulli independientes (v.a. respuesta Y), a cada una de las cuales corresponde una determinada combinación de niveles (x0,x1,…xR) de las R variables explicativas X1,…,XR.

Notación:
  • xq=(xq0,xq1,…,xqR)’  (q=1,…,Q) la q-ésima ciombinación de valores de las R variables explicativas en la muestra.
  • nq es el número de observaciones muestrales (con X=xq), con sum(nq)=N.
  • yq es el número de respuestas Y=1, e Yq es el número de respuestas Y=1 en cada xq.
  • Q es la muestra de v.a. independientes Yq con distribuciones B(nq,pq), donde pq=P[Y=1|X=xq] y E[Yq]=nq*pq.
Tenemos dos casos:
         i.  Q=N: cada individuo muestral tiene una combinación diferente de niveles de las R variables explicativas (1 observacxión de la v.a. respuesta Y en cada combinación).
       ii.    Q1 observación de la v.a. respuesta Y en cada combinación).

Entonces la fórmula del modelo es:
pq= sum(beta_r*xqr)) / (1+e(sum(beta_r*xqr)))
o de modo equivalente: Lq=ln[pq/(1-pq)] = sum(beta_R*xqr)
(en forma matricial: L=X*beta)


III.            Estimación por máxima verosimilitud
Los estimadores de máxima verosimilitud (MV) son los valores de los parámetros que dan máxima probabilidad (verosimilitud) a los datos observados. Para hallarlos hay que maximizar la función de verosimilitud de los datos respecto de los parámetros del modelo logit:
         i.            Estimación MV iterativa con Newton-Raphson
        ii.            Estimación por mínimos cuadrados ponderados
       iii.            Propiedades de los estimadores MV

 IV.            Inferencia en regresión logística: para extrapolar los resultados muestrales a la población
        a.       Contrastes de bondad de ajuste
            b.      Contrastes sobre los parámetros del modelo
            c.       Intervalos de confianza
  
V.            Validación
            a.       Residuos
            b.      Medidas de influencia
            c.       Métodos gráficos

VI.            Selección del modelo más apropiado

sábado, 12 de diciembre de 2009

0
Curvas y superficies de respuesta

Ajuste de curvas y superficies de respuesta

    Los diseños factoriales son muy útiles para el tamizado de factores, es decir, para identificar los factores más importantes que afectan el desempeño de un proceso (o caracterización del proceso, Montgomery, 2004). Una vez que se ha identificado el sunconjunto adecuado de variables del proceso, generalmente el paso siguiente es la optimización del proceso, que implica encontrar el conjunto de condiciones de operación de las variables del proceso que producen el mejor desempeño del mismo.
   Puede resultar útil ajustar una curva de respuesta a los niveles de un factor cuantitativo para que el investigador cuente con una ecuación que relacione la respuesta con el factor. Esta ecuación podría utilizarse para hacer interpolaciones, es decir, para predecir la respuesta en niveles intermedios entre los factores, respecto de los que se utilizaron realmente en el experimento.
   Cuando al menos dos de los factores son cuantitativos, puede ajustarse una superficie de respuesta para predecir   con varias combinaciones de los factores de diseño. En este sentido cabe destacar la metodología de superficies de respuesta (RSM), como el enfoque de optimización más exitoso y generalizado. En general, se usan métodos de regresión lineal para ajustar estos modelos a los datos experimentales. Además, los efectos de los factores cuantitativos pueden representarse con efectos polinomiales con un solo grado de libertad. De manera similar, es posible hacer la partición de las interacciones de factores cuantitativos en componentes de interacción con un solo grado de libertad. 

   El enfoque usual es utilizar el diseño de experimentos para determinar cuáles variables están influenciando la respuesta de interés.  Una vez que dichas variables son identificadas, se obtiene un estimado aproximado de la superficie de respuesta por medio de modelos factoriales especiales.  Esta superficie de respuesta se usa como guía para variar gradualmente los factores controlables que afectan la respuesta de manera tal que se mejore el valor de la respuesta.  Una vez que el cambio de los factores controlables no origine una mejora predecible en la variable de la respuesta, se puede aplicar un método de experimentación más sofisticado para encontrar la superficie de respuesta operativa final del proceso de interés. 
   Supongamos que un investigador desea encontrar los niveles de la variable  x1 y de  x2 que maximicen el rendimiento Y de un proceso. El rendimiento del proceso es una función de los niveles de la variable x1 y x2: Y=f(x1,x2)+epsilon.
   Donde epsilon representa el ruido o error observado en la respuesta Y. Si el valor de la respuesta se denota por E(Y)=f(x1,x2) entonces la superficie representada por  E(Y)=f(x1,x2)se llama superficie de respuesta. Para ayudar a visualizar la forma de una superficie de respuesta, con frecuencia se trazan los contornos de la superficie de respuesta (líneas de respuesta constante en el plano x1,x2). Cada contorno corresponde a una altura particular de la superficie de respuesta. La gráfica de contorno es útil para estudiar los niveles de  x1,x2 que producen cambios en la forma o altura de la superficie de respuesta.
   En la mayoría de los problemas de RMS no se conoce la forma de la relación entre la respuesta y las variables independientes. Por tanto el primer paso es encontrar una aproximación adecuada de la verdadera relación funcional entre   y las variables indepdendientes, por lo general se emplea un polinomio de orden inferior en alguna región de las variables independientes. Si una función de las variables independientes modela adecuadamente la respuesta, entonces a función de aproximación es el modelo lineal de primer orden:
 Y=beta0+beta1*x1+...+betaK*xK+epsilon
   Si hay curvatura en el sistema, entonces se debe utilizar un polinomio de orden superior, tal como el modelo lineal de segundo orden:        
  Y=beta0+sum(beta_j*x_j)+...+sum(beta_ij*x_ij)+epsilon
   El modelo lineal ocasiona que se modele la superficie de respuesta con líneas rectas o planos, mientras que las superficies de respuesta de segundo orden y mayores corresponden a formas geométricas más complejas. En ocasiones se utilizan modelos más complejos para la respuesta de superficie (Montogomery (2000), Montogomery y Myers (1995).

viernes, 11 de diciembre de 2009

0
Modelos lineales generalizados (GLM)

Modelos lineales generalizados

   Los  modelos lineales generalizados (GLM) tienen como objetivo describir el efecto de una o más variables explicativas (independientes) sobres una o más variables respuesta (dependientes).

1 o más variables explicativas    -(efecto)->    1 o más variables respuesta
(X0, …, Xd independientes)                    (Y0, …, Yr dependientes)

   Las componentes del vector Y son variables independientes con distribución proveniente de una familia exponencial. Las variables X0, …, Xd originan un predictor lineal η dado por η =β0 + β1X1 +...+βdXd, o en forma matricial η = Aβ.
   El predictor lineal y la variable dependiente están relacionados por una función de enlace o link g, siendo g monótona y diferenciable. Podemos escribir el modelo en cualquiera de sus 3 formas:
 

  1.   Familia exponencial: una variable aleatoria Y tiene distribución proveniente de una familia exponencial si su función de probabilidad puntual o de densidad es de la forma:



    Si  a(y)=y  se dice que f está dada en su forma canónica. Ejemplos: Poisson, Normal y Binomial



  2. Un GLM queda especificado mediante tres componentes:
    1. Componente aleatoria: distribución de probabilidad de la variable respuesta Y que pertenece a la familia exponencial natural. 
    2. Componente sistemática: función lineal de las variables explicativas que se usa como predictor lineal.
    3. Función de enlace o ligadura: función g que describe la relación funcional entre la componente sistemática y el valor esperado de la componente aleatoria.
               
               Tipos de Modelos Lineales Generalizados para el Análisis Estadístico:

Comp. aleatorio    Link    Comp. sistemático    Modelo
Normal                           Identidad              Continuo                 Regresión
      Normal                           Identidad              Categórico      Análisis de varianza
       Normal                           Identidad              Mixto           Análisis de covarianza
        Binomial                         Logit                    Mixto                 Regresión logística
Poisson                           Log                     Mixto                      Log-lineal
            Multinomial            Logit  generalizado        Mixto                Respuesta multinomial

      3.   Comparación con modelos lineales:





Modelo lineal simple (ML)

   El modelo de probabilidad lineal o modelo de regresión lineal, es de la forma:
 Y(x)=a+βx+epsilon(x)
 donde E[Y|X=x]=p(x)=a+βx. Sin embargo, este modelo no puede explicar el comportamiento de las probabilidades de respuesta de una variable aleatoria binaria:
  • porque las probabilidades pertenecen al intervalo [0,1] y la función lineal toma valores en toda la recta real
  • porque no existe homocedasticidad: la varianza de la variable respuesta Var[Y|X=x]=p(x)(1-p(x)), es decir, no es constante para todos los valores de X (los estimadores mínimos cuatrados son ineficientes)
  • La variable Y no es Normal (no se pueden utilizar distribuciones muestrales de lso estimadores mínimos cuadrados para su inferencia)
  • la tasa de cambio es constante: variaciones iguales de la probabilidad de respuesta frente a variaciones iguales de la variable explicativa. Sin embargo, esto no se cumple para variables respuesta aleatorias binarias.


Modelos no-lineales


  Los modelos no lineales son de la forma:
Y(x)=F(a+βx)+epsilon(x)
donde p(x)=F(a+βx), con F una función distribución estrictamente creciente y F-1(p(x))=a+βx. La relación entre x y p(x) es curvilínea, monótona y acotada entre [0,1].

Ejemplos de F: función de distribución logística para el Modelo de rergesión logística;(1), función de distribución normal para el modelo probit (2) y función de distribución Gumbel para el modelo de valores extremos (3).

  1. Modelo de regresión logística (logit)

   El modelo de regresión logística o función logística, es un modelo lineal generalizado en el cual la variable de respuesta Y es binaria, siendo p la probabilidad de que ocurra el evento en cuestión.
  El modelo simple es de la forma:
p(x)=e(a+βx)/(1+e(a+βx))=1/(1+e-(a+βx))

o lo que es lo mismo, su transformación logit: ln(p(x)/(1-p(x)))=a+βx
   La función de enlace o link adecuada es la función logit (ec.1),  cuya función inversa es la función logística   (ec.2). Si (ec.3), el  modelo de regresión logística está dado por (ec.4), o lo que es lo mismo (ec.5).


Interpretación de los coeficientes:
    La ventaja de la respuesta Y=1 para el calor observado de X=x viene dado por el cociente Odds=p(x)/(1-p(x)).  El Odds es una medida de riesgo que nos indica  cuanto más probable es que ocurra un evento respecto a que no ocurra β0 -> e(β0).
    Entonces se define el riesgo relativo de respuesta Y=1 para dos valores distintos x1 y x2 de la variable explicativa X: R12=p(x1)/p(x2).
    El cociente de ventajas (Odds ratio, OR) de respuesta Y=1 dados x1 y x2 distintos, es:  OR=[p(x1)/(1-p(x1))]/[p(x2)/(1-p(x2))]. Este permite estimar el incremento de “riesgo” de la variable dependiente  por unidad de variación de las variables independientes.  El Odds ratio (OR) cuantifica la magnitud de la relación entre la respuesta y el cambio en una unidad del factor de interés.

   Además, la relación de los cociente de entajas con el riesgo relativo, viene dada por:
OR=R12x(1-p(x2))/(1-p(x1))
   Haciendo cuentas, tenemos que cuando X aumenta una unidad, el incremento de riesgo es OR=e(β1). Si X aumenta k unidades OR=e(kβ1).



  1.  
  2. Modelo probit
  3. Modelo de valores extremos






domingo, 13 de septiembre de 2009

0
aplicaciones

0
jerarquicos

0
Diseños factoriales


    Los diseños factoriales surgen cuando se desea investigar el efecto de dos o más factores sobre una variable respuesta estudiando todas las posibles combinaciones de sus niveles. Por ejemplo, cuando se tienen dos factores A y B; con a y b niveles, respectivamente, se investigan, en cada ensayo completo o réplica del experimento, todas las posibles combinaciones ab de los tratamientos. 
    Para determinar el efecto de cada factor en la respuesta, se debe determinar el cambio experimentado por la respuesta cuando se produce un cambio en el nivel del factor. Este efecto se refiere como efecto principal puesto que se debe a los factores de interés. En cambio, cuando se analizan los cambios o variabilidad de la respuesta debida a la interacción entre dos o más factores, no se habla de efectos principales sino de efectos de interacción. Cuando la magnitud de dichos efectos es elevada, los efectos debidos a factores principales tienen poco significado práctico. Así, en esta situación sería necesario mantener fijos los niveles de los otros factores, para el análisis de un factor principal en presencia de interacciones significativas.
    Los diseños factoriales son más eficientes que los diseños de un factor. Adicionalmente, este tipo de diseños son necesarios cuando hay interacciones significativas entre los factores para evitar conclusiones engañosas. Los diseños factoriales permiten asimismo estimar los efectos de un factor en diferentes niveles de otros factores, obteniéndose resultados e interpretaciones que son válidas bajo todas las condiciones experimentales.



Material (pincha aquí)

  • Diseño factorial (uni-, bi- y tri-factorial)
  • Ajustes de curvas y superficies de respuesta  

Análisis estadístico: modelo bifactorial





En particular, para k=1 réplicas tenemos:


0
cuadrado latino y grecolatino

0
bloques

  En un diseño aleatorizado por bloques completos se consideran tres fuentes de variabilidad: el factor de tratamientos, el factor de bloques y el error aleatorio. La palabra completo se debe a que en cada bloque se prueban todos los tratamientos, es decir, que los bloques están completos. La aleatorización se hace dentro de cada bloque; no se realiza de manera total como en el diseño completamente al azar. 
  Además del diseño completamente aleatorizado por bloques, existen otros diseños interesantes que explotan la idea de la formación de bloques, por ejemplo: el diseño de bloques incompletos balanceados. 
  
  El análisis por bloques permite controlar sistemáticamente fuentes de variabilidad externas que intervienen o actuan en la observación de la variable respuesta. Puesto que en estos casos, la componente aleatoria de error o error experimental asociada a cada observación refleja tanto el error aleatorio como la variabilidad debida a los elementos físicos que intervienen en la realización del experimento. Por tanto, se considerara no sólo el diseño de factores internos sino también la actuación de las fuentes de variabilidad externas que se reflejan en el diseño de los bloques o variación entre bloques.

  • Diseño por bloques aleatorizados completos
    • Adecuación del modelo
    • Estimación mínimo cuadrática de los parámetros y contrastes de significación
  • Diseño por bloques aleatorizados incompletos
    • balanceado
    • Estimación mínimo cuadrática de los parámetros
    • Información inter-bloque en el diseño balanceado

  •    Contrastes de comparación múltiple en el caso del diseño aleatorizado por bloques completos considerando efectos fijos de los tratamientos y bloques.
Si en un diseño aleatorizado por bloques los tratamientos son fijos y el análisis de varianza indica que existe una diferencia significativa entre las medias de tratamiento, el investigador estará interesado en realizar comparaciones adicionales en grupos de medias de tratamiento, para determinar cuáles son las medias que difieren. Con algunas variantes, cualquier método estudiado en el diseño unifactorial puede ser utilizado para este fin. Para llevar a cabo las comparaciones entre grupos de tratamiento para un diseño aleatorizado por bloques se debe sustituir el número de réplicas o repeticiones (n) por el número de bloques (b) en las fórmulas utilizadas en cada uno de los métodos estudiados en los diseños unifactoriales y además se debe utilizar los grados de libertad del error que están definidos por (a-1)(b-1) para un diseño aleatorizado por bloques. A continuación se presentarán los métodos descritos para el diseño unifactorial, expresando solamente las variantes que se deben incorporar para llevar a cabo la comparación de medias de tratamiento para un diseño aleatorizado por bloques. Las hipótesis a probar, el procedimiento y conclusiones se harán de igual manera que para el diseño unifactorial.
  1.   Comparación de Parejas de Medias de Tratamientos.
    1. Método de la Mínima Diferencia Significativa (LSD). El LSD estará dado de la siguiente manera: LSD= t_alfa/2(a-1)(b-1)*sqrt(2MS_e/b)
    2. Prueba de Intervalos Múltiples de Duncan. El error estándar de cada promedio se calcula de la siguiente forma: S_Yi· = sqrt(MS_e/b) . Para encontrar los intervalos significativos r_alfa(p,f), para p=2,3,...,a, a sigue siendo .... el nivel de significancia y f el número de grados de libertad del error que son (a-1)(b-1). De igual manera para  encontrar los mínimos intervalos significativos R_p = r_alfa(p,f)*S_Y·  con p=2,3,...,a, se tomará f como el número de grados de libertad del error (a-1)(b-1).
    3. Prueba de Tukey. El valor crítico de todas las comparaciones vendrá dado por: T_alfa = q_alfa (a,f)*S_Y· donde S_Y· es el error estándar de cada promedio y está dado por S_Y·=sqrt(MS_e/b) y f=(a-1)(b-1) los grados de libertad del error, alfa el nivel de significancia y a el número de tratamientos.
  2.  Comparación de Medias de Tratamientos Individuales 
  3.  Comparación de Tratamientos con un control.
  •    Eficiencia relativa del diseño aleatorizado por bloques. Relación existente entre esta fórmula y la estimación de derivada.
   El análisis o estudio de un experimento puede ser llevado a cabo a través de un diseño unifactorial (o diseño Completamente Aleatorizado) o por un diseño completamente aleatorizado por Bloques. Sin embargo, se puede dar el caso que no se obtenga la misma sensibilidad. En general al utilizar un diseño unifactorial la suma de cuadrados medios del error (MS_e) podría ser mayor que al utilizar un diseño aleatorizado por bloques; ya que el diseño aleatorizado por bloques reduce suficientemente la cantidad de ruido para lograr detectar diferencias significativas entre los tratamientos.
   En un diseño aleatorizado por bloques resulta útil estimar la eficiencia relativa, para compararlo con el diseño unifactorial. El valor que resulta de esta estimación se puede interpretar como el incremento del número de réplicas necesarias que hay que llevar a cabo en un diseño unifactorial para que pueda ser usado en lugar de un diseño aleatorizado por bloques, y así mantener la misma sensibilidad en ambos diseños.

   La forma de definir la eficiencia relativa es mediante la siguiente fórmula:
R=[(df_b+1)*(df_CS+3)*sigma2_CA ]/ [(df_b+3)*(df_CA+1)*sigma2_b]

donde sigma2_CA y sigma2_b son la varianza del error experimental del diseño unifactorial  y del diseño aleatorizado por bloques, respectivamente, y df_CA y df_b sus grados de libertad (df_CA=N-a  y df_b=(a-1)*(b-1)).  Se realiza entonces un ajuste sobre la diferencia grados de libertad entre los dos diseños mediante la razón de grados de libertad en R.
   Como puede observarse para calcular la eficiencia relativa, se deben llevar a cabo estimaciones para sigma2_CA y sigma2_b, las cuales es posible estimarlas de la siguiente forma: sigma2_b ~ MS_e del diseño aleatorizado por bloques, y sigma2_CA=[ (b-1)*MS_bloques + b*(a-1)*MS_e ]/(ab-1), es un estimador insesgado de la varianza del error de un diseño unifactorial, con MS_bloques la suma de cuadrados medios del efecto de los bloques.

  •    Describir cómo es la formulación iterativa de la estimación de un dato faltante en el caso de dos datos faltantes.
   Algunas de las observaciones en uno de los bloques puede hacer falta, cuando se utiliza un Diseño Aleatorizado por Bloques Completos; esto puede suceder debido a algún descuido o error, o por razones fuera de control del experimentador, como la pérdida de alguna unidad experimental. Una observación faltante genera un problema en su análisis, ya que hace que el Diseño este desbalanceado, y se dice que los tratamientos y los bloques no son ortogonales, porque todos los tratamientos no ocurren en todos los bloques. Existen varias formas de solucionar este problema, una de ellas es realizar un análisis aproximado en el que se estima la observación faltante y luego se lleva a cabo el Análisis de Varianza tomando la observación estimada como si fuera un dato real. Este análisis aproximado consiste en hacer estimaciones de los valores faltantes, de manera que se minimice la media de cuadrados del error.
   Supóngase que falta la observación Y_ij que corresponde al tratamiento i y al bloque j; y se representa por x.
   El procedimiento que se lleva a cabo para estimar Y_ij, es el siguiente:
  1. Se calcula el gran total con la observación faltante y se representa por Y`_·· .
  2. Se obtienen los totales del tratamiento y del bloque con el dato faltante que se representa por Y`_i· y Y`_·j respectivamente.
  3. Se calcula el estimador de la observación faltante de la siguiente forma: Y_ij = [ aY`_i· + bY`_·j - Y_·· ] /[(a-1)(b-1)]
   Para más detalles (Ver Douglas C. Montgomery, 1991, Pág 133 y 134)

   Puede suceder que falte más de una observación en el experimento. Existen dos formas para encontrar estas observaciones:
  1. Utilizar el procedimiento descrito anteriormente iterativamente. Por ejemplo, supóngase que hacen falta dos observaciones, la forma de llevar a cabo la estimación de las dos observaciones, es estimando arbitrariamente el primer valor faltante y se usa este valor como un dato real para estimar el segundo. Luego se hace una segunda estimación para el primer dato faltante utilizando la estimación del segundo; con la estimación encontrada para el primero se vuelve a estimar el segundo. Este procedimiento continúa hasta obtener la convergencia en los valores estimados; es decir, hasta que resulten valores parecidos en cada iteración.
  2. Escribir la suma de cuadrados del error en función de los datos faltes, derivar con respecto a cada uno, igualar a cero y resolver las ecuaciones que resultan. 
   En general, para cualquier problema que falten datos, el número de los grados de libertad del error se debe reducir en uno por cada dato que es estimado.

   Además del diseño completamente aleatorizado por bloques, existen otros diseños interesantes que se introducen brevemente a continuación y que explotan la idea de la formación de bloques, por ejemplo: el diseño de bloques incompletos balanceados.


  •   Diseños por bloques incompletos balanceados
   Puede darse el caso que en algunos experimentos que usan Diseños Aleatorizados por Bloques no se puedan llevar a cabo los ensayos de todas las combinaciones de tratamientos dentro de cada bloque; ya sea por la escasez de los recursos del experimento, por la situación económica, o por el tamaño físico de los bloques. Para analizar estos tipos de experimentos se usa el Diseño Aleatorizado por Bloques en los que cada tratamiento no está presente en cada bloque; y a este tipo de Diseño se conocen como Diseño Aleatorizado por Bloques Incompletos.
   Un Diseño particular de ellos son los Diseños por Bloques Incompletos Balanceados; el cual consiste en un Diseño por Bloques Incompleto en el que cualquier par de tratamientos ocurren juntos el mismo número de veces. Si se tienen tratamientos a y se pueden probar k (k tratamientos en cada bloque entonces un Diseño Balanceado por Bloques Incompletos puede ser construido tomando (a k) combinaciones de bloques y asignándose una combinación de tratamientos diferentes a cada bloque. Sin embargo frecuentemente es posible obtener un Diseño Balanceado con menos de (a k) combinaciones de bloques.
  • REPRESENTACIÓN SIMBÓLICA DE LOS DATOS
  •  MODELO ESTADÍSTICA
  • SUMAS Y MEDIAS DE CUADRADOS
  • ANÁLISIS ESTADÍSTICA
  • ESTIMACIÓN DE LOS PARÁMETROS DEL MODELO
  • COMPARACIÓN ENTRE TRATAMIENTOS
  •   Análisis de datos-Ejemplo
DISEÑO COMPLETAMENTE ALEATORIZADO PRO BLOQUES en R
 
> Block=rep(c("P15","P20","P25","P30","P35"),c(4,4,4,4,4))
> Treat=rep(c("A","B","C","D"),5)
> Datos=c(7,12,14,19,7,17,18,25,15,12,18,22,11,18,19,19,18,19,23,11)
> Data=cbind(Datos,Treat,Block)
> par(mfrow=c(1,2),cex=.8)
> interaction.plot(Treat,Block,Datos,type="b",legend=FALSE)
> interaction.plot(Block,Treat,Datos,type="b",legend=FALSE)
> M<-matrix(Datos,nrow=4,ncol=5)
> dimnames(M)=NULL
> rownames(M)=c("A","B","C","D");colnames(M)=c("P15","P20","P25","P30","P35")
> M
  P15 P20 P25 P30 P35
A   7   7  15  11  18
B  12  17  12  18  19
C  14  18  18  19  23
D  19  25  22  19  11

> library(lattice)
> B=stripplot(Block~Datos|Treat,layout=c(4,1))
> A=stripplot(Treat~Datos|Block,layout=c(5,1))
> print(A,plit=c(1,1,1,2),more=T)
> print(B,plit=c(1,2,1,2),more=F)

> Data=as.data.frame(Data)
> Da=as.list(Data)
> plot.design(Datos~Da$Treat+Da$Block)
> mod.aov=aov(Datos~Da$Treat+Da$Block)
> summary(mod.aov)

> library(PASWR)
> checking.plots(mod.aov)
> CI=TukeyHSD(mod.aov,which="Da$Treat")
> Media=NULL;for(j in 1:4){Med=mean(M[j,]);Media=c(Media,Med)}
> barplot(Media,ylim=c(0,30))
> library(psych)
> error.bars(t(M),add=T)
> par(mfrow=c(2,2))
> plot(mod.aov,pch=20)

> Mm=NULL;for(i in 1:4){S=mean(M[i,]);Mm=c(Mm,S)}
> Mn=NULL;for(i in 1:5){S=mean(M[,i]);Mn=c(Mn,S)}
> m=cbind(M,Mm)
> m=rbind(m,c(Mn,mean(Mm)))
> SSblock=4*sum((m[5,1:5]-mean(Datos))^2)
> SStrat=5*sum((m[1:4,6]-mean(Datos))^2)
> SS=NULL;for(j in 1:4){S=sum(((M[j,]-mean(M)))^2);SS=c(SS,S)}
> SStot=sum(SS)
> SSerror=SStot-SStrat-SSblock

    0
    regresion lineal


    smooth relationships (Bolker 2007_ Ecological Models and Data in R)

    • ˆ R incorporates two slightly different versions of robust locally weightedregression (lowess and loess). This algorithm runs linear or quadratic regressions on successive chunks of the data to produce a smooth curve. lowess has an adjustable smoothness parameter (in this case the proportion of points included in the “neighborhood” of each point when smoothing) that lets you choose curves ranging from smooth lines that ignore a lot of the variation in the data to wiggly lines that pass through every point: in Figure 2.8a, I used the default value (lines(lowess(Initial,Killed))).
    • ˆ Figure 2.8a also shows a spline fit to the data which uses a series of cubic curves to fit the data. Splines also have a smoothing parameter, the degrees of freedom or number of different piecewise curves fitted to the data; in this case I set the degrees of freedom to 5 (the default here would be 2) to get a slightly more wiggly curve (smooth.spline(Initial, Killed,df = 5)).
    • ˆ Simpler possibilities include just drawing a straight line between the mean values for each initial density (using tapply(Killed,Initial,mean) to calculate the means and unique(Initial) to get the non-repeated values of the initial density), or plotting the results of a linear or quadratic regression of the data (not shown: see the R supplement). I plotted straight lines between the means in Figure 2.8b because local robust regression and splines worked poorly.


    1
    Análisis de varianza: ANOVA. Efectos fijos y aleatorios en R.

    En esta entrada evaluaremos el análisis de varianza con efectos fijos y efectos aleatorios, así como aquellas situaciones donde tenemos medidas con repeticiones. Como resumen teórico:
    Resumen Anova Efectos Fijos y Aleatorios


    Nuestros objetivos son:
    • Deducir los valores esperados de las sumas de cuadrados medios (MS) entre tratamientos y dentro de cada tratamiento para el modelo de efectos fijos y efectos aleatorios.
    • Evaluar cuál de los cuatro métodos estudiados para realizar comparaciones múltiples entre niveles medios de tratamientos por pares posee mejores propiedades.
    • ¿Cómo afecta la violación de la hipótesis de normalidad a la prueba F en los modelos de efectos fijos y aleatorios?. ¿Influye de forma distinta en diseños con tamaños muestrales fijos y variables?.
    • Elaborar un resumen sobre los métodos no paramétricos usuales en el análisis de varianza.
    Aquí les dejo un resumen teórico del tema completo, y un ejemplo aplicado en R
    anova_1_via

    • Ejemplo: Análisis de datosen R (r-project)
    ########### ejemplo: ANOVA ##############
    ##Data
    A_15=c(7,7,15,11,9)
    B_20=c(12,17,12,18,18)
    C_25=c(14,18,18,19,19)
    D_30=c(19,25,22,19,23)
    E_35=c(7,10,11,15,11)
    ##Examen gráfico
    scores=data.frame(A_15,B_20,C_25,D_30,E_35)
    boxplot(scores)
    library(PASWR)
    scores2=stack(scores) #preparación de los datos
    X<-scores2[,1]
    INDEX<-scores2[,2]
    oneway.plots(X,INDEX) #dotplot, boxplot y
    design plot (means)

    ## Modelo con efectos fijos (FIXED MODEL)
    # Las medias de los tratamientos son o no iguales: Ho: mu1=mu2=...=mua vs Ha: mui!=muj
    # Cuando la Ho es cierta, se puede evaluar una afirmación equivalente en término de los efectos
    de los tratamientos: Ho: tau1=tau2=...=taua vs Ha: taui!=0
    ##Estimaciones: E(MSerror)=sigma^2 y E(MStrat)=sigma^2 + sum((ni*taui^2)/(a-1))
    #Prueba:
    TreatMean<-tapply(X,INDEX,FUN=mean)
    a<-length(TreatmentMean)
    N<-length(X)
    dft<-a-1
    dfe<-N-a
    GrandMean<-mean(X)
    ni<-nrow(scores)
    SStreat<-ni*sum((TreatMean-GrandMean)^2)
    SStot<-sum((X-GrandMean)^2)
    SSerror<-SStot-SStreat
    MStreat<-SStreat/dft
    MSerror<-SSerror/dfe
    Fobs<-MStreat/MSerror
    pvalue<-1-pf(Fobs,dft,dfe)
    ##equivalente
    summary(aov(X~INDEX))
    model.tables(aov(X~INDEX),type="means")
    ##Chequeo de supuestos; 3 supuestos en el componente de ERROR (se util los residuales como su estimador): independencia, distribución normal y varianza constante
    mod.aov<-aov(X~INDEX)
    library(MASS)
    r<-stdres(mod.aov)
    n<-length(X)
    #Chequeo de independencia de errores
    par(pty="s")
    plot(1:n,r,ylab="Standardized
    Residual"
    ,xlab="Ordered Value")
    #Chequeo de normalidad: qqplot y shapiro.test()
    par(pty="s")
    qqnorm(r)
    abline(a=0,b=1)
    shapiro.test(r)
    #Chequeo de varianza constante: plot de residuos (estandarizdos) vs valores ajustados; Levene
    tm<-fitted(mod.aov)
    plot(tm,r,xlab="Fitted Value",ylab="Standardized Residual")
    med<-tapply(X,INDEX,median)
    ZIJ<-abs(X-med[INDEX])
    summary(aov(ZIJ~INDEX))
    library(PASWR)
    checking.plots(mod.aov) ##package PASWR
    ##comparación múltiple de medias
    #Ho: Ho1 intersección Ho2 intersección ...HoK
    library(multcomp)
    library(multcompView)
    CI<-TukeyHSD(aov(X~INDEX,which="INDEX"))
    plot(CI,las=1)
    INDEX.aov<-aov(X~INDEX)
    MSE<-summary(aov(INDEX.aov))[[1]][2,3]
    alpha.c<-0.05
    ybari<-TreatmentMean
    TcritLSD<-qt(1-alpha.c/2,dfe)
    nn<-rep(ni,a)
    LSD<-TcritLSD*sqrt(MSE)*sqrt(sum(1/nn))
    TcritTUK<-qtukey(1-alpha.c/2,a,dfe)/sqrt(2)
    HSD<-TcritTUK*sqrt(MSE)*sqrt(sum(1/nn)) #nn es un vector de ni y nj, con el length=número de tratamientos
    library(gregmisc)
    NS<-tapply(X,INDEX,length)
    SE<-sqrt(MSE)/sqrt(NS)
    t.v<-qt(.95,dfe)
    ci.l<-ybari-t.v*SE
    ci.u<-ybari+t.v*SE
    barplot2(ybari,plot.ci=T,ci.l=ci.l,ci.u=ci.u,col="sky blue",ci.lwd=2)
    title(main="Mean X por INDEX \n con CI individual 95%")
    #multcompBoxplot(X~INDEX) Su gráfico no es fácilmente interpretable

    ## Modelo con efectos aleatorios (RANDOM MODEL)
    ##supuestos: eijNID(0,sigma), taui~NID(0,sigma), taui y eij son independientes
    #Ho: sigma-subtau^2=0 vs Ha: sigmasubtau^2>0
    #Estimaciones:
    #cuando los a tratamientos tienen igual tam de muestreo: sig2=estim(sigma)^2=MSerror y sig2tau=estim(sigma)-subtau^2=(MStreat- MSerror)/n
    #cuando los tamaños muestrales son desiguales, n se reemplaza por n`=1/(a-1)*sum(ni)-(sum(ni^2)/sum(ni))
    summary(aov(X~INDEX))
    MSC<-summary(aov(X~INDEX))[[1]][1,3]
    MSE<-summary(aov(X~INDEX))[[1]][2,3]
    #Estimación de los componentes de varianza
    sig2tau<-(MSC-MSE)/n #nº de tratamientos