---
title: "Modelo lineal normal"
author: "Guillermo Ayala Gallego"  
date: "`r Sys.Date()`"  
address: Departamento de Estadística e I.O. Universidad de Valencia  
format: 
    revealjs:  
      scrollable: true
      echo: true
      embed-resources: true
    pdf:
       toc: true
       number-sections: false
       colorlinks: true  
---


```{r}
#| label: edm0 
#| echo: false
dirTamiData = "/home/gag/ownCloud/alltami/tami-data/"
dirTamiFigures = "/home/gag/ownCloud/alltami/tami13/manual/figures/"
```


# Introducción  

## ¿Qué información tenemos?

-  Tenemos información (numérica, categórica) sobre una serie
de muestras u observaciones.
- Para la $i$-ésima observación tenemos $p$ variables que recogemos en
el vector $\mathbf x_i \in \mathbb R^p$ donde
$$
\mathbf x_i =  
\begin{bmatrix}
x_{i1}\\
\vdots\\
x_{ip}
\end{bmatrix}.
$$
-  $\mathbf x_i^T = (x_{i1},\ldots,x_{ip})$.
-  Estas variables pueden ser númericas o categóricas.
-  Pueden valores fijados por nosotros (dosis de medicación por ejemplo)
o bien observados y por lo tanto realizaciones de variables aleatorias.

## Variable respuesta  

- En lo que sigue   consideramos que hay una variable **importante**.     
- La variable importante (para el experimentador) se le llama **variable respuesta**.
- En denominaciones más clásicas, **variable dependiente**.  
- El resto de variables del estudio (las que recogemos en $\mathbf x_i$) son las
variables predictoras o independientes.  
- Nuestra información está formada por $(\mathbf x_i, y_i)$ con $i=1,\ldots,n$.  

## ¿Qué problema queremos resolver?
	
- Una respuesta fácil es decir que queremos conocer el valor de
  la variable respuesta utilizando las variables predictoras.
  - La respuesta anterior es **falsamente** simple.
  - ¿Solamente queremos conocer el valor de la respuesta?
  - Quizás estamos pensando en un futuro en donde conozcamos las variables
  predictoras  y nos interese saber cuál será la respuesta correspondiente.
  - ¿El valor exacto de la respuesta?
  - ¿La media de la respuesta?
  - ¿Un valor numérico que aproxime cada una de estas cantidades o bien un
	intervalo que las contenga?

## Media condicionada
	
- Los valores observados $y_i$ consideraremos que son realizaciones
  de una variable aleatoria $Y$.
  - Realmente en lo que sigue modelizaremos el comportamiento aleatorio
  de la variable $Y_i$ condicionada a los valores observados $\mathbf x_i$.
  - Es decir, nuestro interés estará en **la distribución   condicionada**.
  -  Los valores $Y_i$ serán independientes entre sí.
  - El interés se centrará (fundamentalmente) en la dependencia de la
  media de $Y_i$ respecto de las covariables.
  - Nuestro interés
  fundamental (pero no único) estará en conocer las  medias condicionadas 
	  $\mu_i = E[Y_i|\mathbf x_i]$. 

## Datos de Galton
	
- Históricamente podemos considerar el origen de este tipo de modelos.
- El problema que se plantea Galton es estudiar la posible relación que 
  puede tener la estatura de un hijo o hija en edad adulta con las estaturas
  de sus padres.

## Datos de Galton$\ldots$  

::: {layout-ncol=2}

![](../figures/GaltonFamiliesScatterplot.png)  

![](../figures/Galton_notebook_page_1.jpg)  

:::

## Modelos sobre la media

-  Nos planteamos conocer la media de la respuesta aleatoria.
- Pretendemos conocer la media de la respuesta aleatoria
**condicionada**  a los valores de las variables predictoras.
- Denotamos
$$
\mu_i = E[Y_i | \mathbf x_i].
$$

# Dependencia lineal  

## ¿Qué tipos de dependencia vamos a considerar?  

-  En la mayor parte de los casos dependencias de tipo lineal.
- Suponemos dos predictores $\mathbf x = (x_1,x_2)^T$ tales que $x_1$ es 
numérico y el segundo es una variable categórica binaria codificada
con 1 y 0.
- ¿Cómo modelizamos la dependencia de la media condicionada respecto de
$x_1$?
- Una dependencia lineal vendría dada como
$$
  \mu_i = \beta_0 + \beta_1 x_{i1}. 
$$

## 

![](../figures/fig311218a.png)


## ¿Y la dependencia de las $\mu_i$ respecto de la variable binaria?

- Obviamente simplemente tenemos dos valores.
- Un modo simple es 
$$
\mu_i = \beta_0 + \beta_2 x_{i2}. 
$$
- Cuando $x_{i2}=0$ 
  $$\mu_i =\beta_0.$$
- Cuando $x_{i2}=1$ 
  $$\mu_i =\beta_0 + \beta_2.$$



## ¿Y las dos variables predictoras conjuntamente consideradas?  

- El modelo más sencillo sería:
$$
\mu_i = \beta_0 + \beta_1 x_{i1} +   \beta_2 x_{i2}.
$$
-  Cuando $x_{i2} = 0$ entonces
$$\mu_i =\beta_0  + \beta_1 x_{i1}. $$
-  Cuando $x_{i2} = 1$ entonces
$$\mu_i =\beta_0 + \beta_2  + \beta_1 x_{i1}.$$

##  
![](../figures/fig311218b.png)


## ¿Como podemos expresar la dependencia de ambas covariables? 

- Otra vez recurrimos a la opción de expresarlo de un modo lineal.
$$
\mu_i = \beta_0 + \beta_1 x_{i1} +   \beta_2 x_{i2} + \beta_3 x_{i1}x_{i2}
$$
- Cuando $x_{i2} = 0$ tenemos 
$$
\mu_i = \beta_0 + \beta_1 x_{i1}. 
$$
Cuando $x_{i2} = 1$ tenemos 
$$
\mu_i = (\beta_0 + \beta_2) + (\beta_1 + \beta_3) x_{i1}. 
$$

##   

![](../figures/fig311218f.png)


## Predictor categórico y variables dummy

- Cuando consideramos una variable predictora categórica con más de dos 
categorías entonces es habitual codificarla utilizando variables
tontas.
- Si la variable predictora categórica tiene $I$ categorías entonces 
se elige una categoría de referencia (por ejemplo, la primera).
-  Las variables binarias asociadas a la variable original $x$
serían: $v_1 =1$ si $x=2$ y cero en otro caso; $v_2 =1$ si $x=3$ y cero
en otro caso; $\ldots$; $v_{I-1} =1$ si $x=I$ y cero en otro caso.
- Obviamente cuando todas las variables
$v$ son nulas estamos en la primera categoría.

##  

- Si solamente tenemos la 
variable categórica como predictora entonces la media sería función lineal
de $x$ del siguiente modo:
$$
\mu_i = \beta_0 + \beta_1 v_1 + \ldots + \beta_{I-1} v_{I-1}.
$$

##  
   
- Por ejemplo, supongamos una numérica $x$ y una categórica $v$ con
$I$ categorías. Construimos las variables tontas siendo $I$ la de referencia.
Podemos considerar modelos como
$$
\mu_i = \beta_0 + \beta_1 x_{i},
$$
que lo expresamos como
$$
y \sim x
$$
- Un modelo que contiene solamente a $v$ sería el dado previamente
$$
\mu_i = \beta_0 + \beta_1 v_{i1} + \ldots + \beta_{I-1} v_{i,I-1}.
$$
y lo expresamos como
$$
y \sim v
$$

##  

   
- Un modelo que contempla ambas variables puede ser 
$$
  \mu_i = \beta_0 + \beta_1 x_{i} +
   \beta_2 v_{i1} + \ldots + \beta_{I} v_{i,I-1}
$$
Este modelo lo podemos abreviar como
$$
y \sim x + v
$$
- Un modelo más completo que contempla la posible interacción sería
\begin{multline*}
  \mu_i =  \beta_0 + \beta_1 x_{i} + 
  \beta_2 v_{i1} + \ldots + \beta_{I} v_{i,I-1} + \\
  \beta_{I+1}  x_{i}v_{i1} + \ldots + \beta_{2I-1}  x_{i}v_{i,(I-1)} 
\end{multline*}
Esto lo indicaremos como
$$
y \sim x * v
$$
$$
y \sim x + v + x:v
$$


# Recta de mínimos cuadrados

##  
  
-  Nuestros datos son $(x_i,y_i)$ con $i=1,\ldots,n$.
- Un modelo sencillo podría ser

$$
    y_i  =  \beta_0 + \beta_1 x_i, \ i = 1, \ldots,n
$$

- No es posible.
- Una  **buena**  solución aproximada.
- Nos conformamos con
  
$$
y_i  \approx   \beta_0 + \beta_1 x_i, \ i = 1, \ldots,n
$$


##  
   
- Una posibilidad para determinar unos **buenos** valores para
     $\beta_0$ y $\beta_1$ es considerar

$$
S_t = \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)^2.
$$

- Podemos obtener los valores que minimizan la función anterior
igualando las derivadas parciales a cero. 

$$
  \frac{\partial S_t}{\partial \beta_0} = -2 \sum_{i=1}^n(y_i -
  \beta_0 - \beta_1 x_i),
$$

$$
  \frac{\partial S_t}{\partial \beta_1} = -2 \sum_{i=1}^n x_i (y_i -
  \beta_0 - \beta_1 x_i),
  $$
   
   
##  
   
- Igualando a cero
$$
  \sum_{i=1}^n y_i - n \hat{\beta}_0 - \hat{\beta}_1  \sum_{i=1}^n x_i = 0,
$$
$$
  \sum_{i=1}^n x_i y_i -  \hat{\beta}_0 \sum_{i=1}^n x_i - \hat{\beta}_1
  \sum_{i=1}^n x_i^2 = 0,
$$

##   
  
- Haciendo algún cálculo se obtiene
  
$$
\hat{\beta}_1 = \frac{\sum_{i=1}^n x_i y_i - (\sum_{i=1}^nx_i \sum_{i=1}^n
y_i)/n}{\sum_{i=1}^n x_i^2 - (\sum_{i=1}^nx_i)^2/n}.
$$

- Se sigue que
	
 $$
\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}.
$$

## Recta de mínimos cuadrados $\ldots$  
  
  - Si denotamos
		
  $$
	  S_{xy} = \sum_{i=1}^n (x_i-\bar{x}) (y_i - \bar{y}),
 $$

y 
		
$$
S_{xx} = \sum_{i=1}^n (x_i-\bar{x})^2, \ \ 
S_{yy} = \sum_{i=1}^n (y_i-\bar{y})^2,
$$
		
entonces 

$$
\hat{\beta}_1 = \frac{S_{xy}}{S_{xx}}.
$$
			
## Predicción y residuo  
  
-  Nos planteamos encontrar unas buenas aproximaciones para 
    $y_i$ utilizando   $x_i$.
- Lo natural sería tomar el siguiente valor la  **predicción**
    o  **valor ajustado** para $y_i$

$$
    \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i.
$$
- El  **residuo** es 

$$
e_i = y_i - \hat{y}_i.
$$ 

# Modelo  

## Regresión lineal simple  

- $Y_i = \beta_0 + \beta_1 x_i + \epsilon_i$ para $i=1,\ldots,n$.
- $\epsilon_i \sim N(0,\sigma^2)$.
- Los errores $\epsilon_i$ son independientes. 

## Regresión lineal simple$\ldots$  
  
- Estamos asumiendo que la distribución **condicionada**
    de $Y_i$ al predictor $x_i$ es normal con media $\beta_0 + \beta_1 x_i$
    y varianza $\sigma^2$, 

$$
Y_i | x_i \sim N(\beta_0 + \beta_1 x_i,\sigma^2),
$$

y que los distintos $Y_i$ son independientes entre sí. 
 

## Representación matricial  

   
- Tenemos  el **vector de respuestas aleatorias**,
     $\mathbf Y$; la **matriz de modelo**,
     $\mathbf X$; el **vector de coeficientes**, $\mathbf \beta$ y el 
     **vector de errores aleatorios**, $\mathbf \epsilon$ dados por
	 
$$
\mathbf Y = 
\begin{bmatrix}
   Y_1 \\
   \vdots \\
    Y_n
\end{bmatrix};
\ \ 
\mathbf X = 
\begin{bmatrix}
  1 & x_1 \\
  \vdots&  \vdots \\
    1 & x_n
\end{bmatrix};
\ \ 
\mathbf \beta = 
\begin{bmatrix}
 \beta_0 \\
  \beta_1
\end{bmatrix};
\ \ 
\mathbf \epsilon = 
\begin{bmatrix}
   \epsilon_1 \\
   \vdots \\
    \epsilon_n
\end{bmatrix}
$$
- El modelo de regresión lineal simple se puede formular como
$$
  \mathbf Y = \mathbf X \mathbf \beta + \mathbf \epsilon.
  $$
   - El vector de errores aleatorios $\mathbf \epsilon$ sigue una
     distribución normal multivariante 
$$
  \mathbf \epsilon \sim N_2(\mathbf 0_p,\sigma^2 \mathbf I_p).
$$
   


## Verosimilitud  
  
  - La función de verosimilitud viene dada por
    \begin{multline} 
      L(\beta_0,\beta_1,\sigma^2;y_1,\ldots,y_n) =
      \prod_{i=1}^n \frac{1}{\sqrt{2\pi} \sigma}
      e^{-\frac{1}{2 \sigma^2} (y_i - \beta_0 - \beta_1 x_i)^2} = \\
      \frac{1}{(2 \pi)^{n/2} \sigma^n}
      e^{-  \frac{1}{2 \sigma^2} \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)^2} 
    \end{multline}
  -  Los estimadores máximo verosímiles de los coeficientes corresponden
    con los estimadores mínimo cuadráticos.
  

# Ejemplo  

## Paquetes  

- Cargamos los paquetes que necesitamos.  

```{r}
#| label: rls1
pacman::p_load(Biobase,ggplot2)
```

## Datos  
- Leemos los datos.  

- La primera es la más ortodoxa. 
```{r}
#| label: rls2
data(gse25171,package="tamidata2")
sel0 = which("261892_at"==fData(gse25171)[,"PROBEID"])
df0 = data.frame(pData(gse25171)[,c("time","Pi")],
                 expression=exprs(gse25171)[sel0,])
```
```{r}
#| label: rls4
#| fig-cap: Expression vs time
#| output-location: slide
(p = ggplot(df0,aes(x=time,y=expression)) + geom_point())
```
- Podemos querer guardar el gráfico.   

```{r}
#| label: rls5
#| eval: false  
ggsave("gse25171_261892_at.png",p)
```
## Recta de mínimos cuadrados  

- Ajustamos. 

```{r}
#| label: rls6
fit = lm(expression ~ time, data = df0)
```  

- ¿Qué clase de objeto tenemos?
```{r}
#| label: rls7
class(fit)
```
-  Vemos lo que el método `print` nos muestra de un objeto de clase
`lm`.
```{r}
#| label: rls8
fit
```
- Los coeficientes de la recta de mínimos cuadrados los podemos obtener
con el método  `coef` o  `coefficients`.

```{r}
#| label: rls9
coef(fit)
coefficients(fit)
```
- Añadimos la recta de mínimos cuadrados al diagrama de puntos.  
```{r}
#| label: rls10
(p = ggplot(df0,aes(x=time,y=expression))+geom_point() +
    geom_smooth(method='lm', se = FALSE))

```
- El método `predict` aplicado al objeto de clase `lm` nos proporciona
las predicciones para las observaciones. 

```{r}
#| label: rls11
predict(fit)
```
- El método `residuals` o `resid` nos proporciona los residuos. 

```{r}
#| label: rls12
residuals(fit)
resid(fit)
```
- El método `summary` aplicado a un objeto `lm`.  
```{r}
#| label: rls13
(fit.s = summary(fit))
```
-  Tenemos un objeto de clase

```{r}
#| label: rls14
class(fit.s)
```

- Algo de la información que contiene el objeto es 
```{r}
#| label: rls16
attributes(fit.s)
```
- La matriz de modelo se obtiene con 
```{r}
#| label: rls18
model.matrix(fit)
``` 

- En el `summary` tenemos los contrastes.

```{r}
#| label: rls19
summary(fit)
```  

## Análisis de la varianza de una vía

- Tenemos una variable de interés $Y$ y pretendemos 
estudiar su posible dependencia de un factor experimental.
- El experimentador tiene un **factor** de interés con distintos niveles
y se pretende evaluar su influencia en la variable respuesta.
- Si el factor tiene dos niveles entonces podemos comparar las medias
mediante un test de la t (t-test).
- Con más de dos niveles del factor hemos de utilizar otras herramientas.
- Nos ocupamos (de un modo muy simple) de lo que se conoce
como **experimentos con un solo factor completamente aleatorizado**.

## Comparando grupos

- Tenemos $I$ condiciones distintas.
- En cada una de ellas $n_i$ muestras.
- $\sum_{i=1}^I  n_i = n$ es  el total de muestras.
- $Y_{ij}$ denota la respuesta aleatoria en la $j$-ésima muestra de la
$i$-ésima condición.

## Modelo

- Suponemos que $Y_{ij}$ con $j=1,\ldots,n_i$ son independientes 
y con la misma distribución.
-  El modelo de análisis de la varianza de una vía es 
\begin{equation}
\label{eq:040219b}
Y_{ij} = \beta_0 + \beta_i + \epsilon_{ij},  
\end{equation}
- Donde se asume que $\epsilon_{ij} \sim N(0,\sigma^2)$ y son
independientes entre si para los distintos grupos y dentro de cada grupo.
- Estamos asumiendo que $Y_{ij} \sim N(\beta_0 + \beta_i,\sigma^2)$.

## Interpretación de los parámetros  

- $\beta_0$ es la media global de todos los grupos. 
- $\beta_i$ sería la diferencia de la media del grupo $i$ respecto de esta
media global.
- $\epsilon_{ij}$ sería el error aleatorio de la observación $(i,j)$ 
respecto del valor medio $\beta_0+\beta_i$. 

##  
   - Se trata de evaluar si hay diferencias entre grupos.
   - Bajo el modelo que
acabamos de formular se traduce en la siguiente hipótesis nula 
$$
H_0: \beta_1 = \beta_2 = \ldots = \beta_I = 0.
$$
frente a que algún par  de los $\beta_i$ sean distintos.


## Contraste  

- ¿Cómo podemos contrastar la hipótesis nula anterior?
- Denotamos por $y_{ij}$ la j-ésima muestra observada bajo la
condición $i$ ($i=1,\ldots,I$ y $j=1,\ldots,n_i$).
- Denotamos las medias muestrales para cada grupo como
$$
\bar{y}_{i \cdot} = \sum_{j=1}^{n_i} \frac{y_{ij}}{n_i}.
$$
- La media de todas las observaciones o media total como  
$$
\bar{y}_{\cdot \cdot} = \sum_{i=1}^I \sum_{j=1}^{n_i} \frac{y_{ij}}{n}.
$$


##   

- Definimos la **suma de cuadrados intra** o **del error**  como 
$$
SS(Within) = \sum_{i=1}^I  \sum_{j=1}^{n_i} (y_{ij} - \bar{y}_{i \cdot})^2,  
$$
y la **suma de cuadrados entre**  como 
$$
SS(Between) = \sum_{i=1}^I n_i (\bar{y}_{i \cdot} - \bar{y}_{\cdot \cdot})^2.  
$$
- 
El estadístico para contrastar esta hipótesis nula es 
$$
F = \frac{SS(Between) / (I-1)}{SS(Within)/(n-I)}.
$$


## Tabla de análisis de la varianza  
| Source | SS | df | MS | F | p |
|--------|----|----|----|---|---|
| Between | SS(B)|I-1| $\frac{SS(B)}{I-1}$ | $\frac{SS(B)/ (I-1)}{SS(W)/(n-I)}$| $P(> F)$ |
| Within | SS(W)| n-I| $\frac{SS(W)}{n - I}$ | |||
| Total  |       SS(B) + SS(W) | |||   
-------



## Contraste$\ldots$  

-  **Bajo la hipótesis nula de que todas las medias son la misma**
(y puesto que asumimos una misma varianza) tendríamos una distribución 
común bajo todas las condiciones.
-  Asumiendo la hipótesis nula el estadístico
$F$ se distribuye como un $F$ con $I-1$ y $n-I$ grados de libertad,
$$
F \sim F_{I-1,n-I}.
$$
- Bajo la hipótesis alternativa, los valores de $F$ tenderán a
ser **grandes** o mayores que los esperables bajo la hipótesis nula.
- La **región crítica** (donde rechazamos la hipótesis nula)
será un intervalo de la forma $[c, +\infty)$.
- Si tomamos como valor $c$ el valor observado tendremos el p-valor.

## Otra formulación del modelo  

-   Consideramos las variables que indican las categorías y que
consideran como categoría de referencia el primer grupo.
- $$
E[Y_{1j}] = \beta_0 
$$
y 
$$
E[Y_{ij}] = \beta_0 + \beta_i
$$
para $i=2, \ldots, I$. De un modo conjunto:
$$
E[Y_{ij}] = \beta_0 + \beta_2 v_{2j} + \ldots + \beta_I v_{Ij}
$$
donde $v_{ij} = 1$ si estamos en el grupo $i$ y cero en otro caso.

- La hipótesis nula de que no hay diferencias entre
las medias de los distintos grupos vendría formulada como
$$
H_0: \beta_2 = \ldots = \beta_I =0.
$$ 


## tamidata2::gse25171  
- Analizamos los datos correspondientes a la sonda `261892_at`.

- Consideramos las variables fenotípicas `time2` y `Pi`. 


```{r}
#| label: rls21
pacman::p_load(Biobase)
data(gse25171,package="tamidata2")
head(pData(gse25171),n=2)
```


```{r}
#| label: rls22
sel0 = which("261892_at"==fData(gse25171)[,"PROBEID"])
df0 = data.frame(pData(gse25171)[,c("time","Pi")],
                 expression=exprs(gse25171)[sel0,])
```

- Construimos una variable categórica combinando las dos predictoras previas. 

```{r}
#| label: rls23
time2Pi = vector("list",ncol(gse25171))
for(i in seq_along(time2Pi))
  time2Pi[[i]] = paste0(pData(gse25171)[,"time2"][i],
                          pData(gse25171)[,"Pi"][i])
time2Pi = factor(unlist(time2Pi))
levels(time2Pi)
```
- Construimos un `data.frame`  en el que consideramos la expresión
de la sonda y la variable que acabamos de construir. 

```{r}
#| label: rls24
sel0 = which("261892_at"==fData(gse25171)[,"PROBEID"])
df1 = data.frame(time2Pi,expression=exprs(gse25171)[sel0,])
summary(df1[,"time2Pi"])
```
-  Realizamos el ajuste del modelo y vemos la tabla anova. 

```{r}
fit2 = lm(expression ~ time2Pi,data=df1)
summary(fit2)
```

```{r}
#| label: rls25
fit2 = aov(expression ~ time2Pi,data=df1)
summary(fit2)
```


#  Mínimos cuadrados  

##   

- Disponemos de los datos $( \mathbf x_i,y_i)$ con $i=1, \ldots, n$ siendo
$\mathbf y = (y_1,\ldots,y_n)^T$, las respuestas observadas;
$\mathbf x_i = (x_{i1},\ldots,x_{ip})^T$, los predictores correspondientes
a la $i$-ésima observación.
- La matriz de modelo que recoge los valores de los predictores:
$$
 \mathbf{X} = 
\begin{bmatrix}
x_{11} & \ldots & x_{1p} \\
\vdots & \vdots & \vdots \\
x_{n1} & \ldots & x_{np} 
\end{bmatrix}
$$
- Asumimos
$$
 \mathbf{\mu} =  \mathbf X  \mathbf{\beta}. 
$$
- ¿Cómo estimamos $\mathbf{\beta}$?


# Planteamiento

##  Estimadores mínimo cuadráticos
-  Una opción clásica son 
los mínimos cuadrados en donde pretendemos que 
$$ 
\|  \mathbf y - \hat{ \mathbf \mu} \|^2
$$
- Lo que estamos haciendo es sustituir el vector
de medias desconocido por los valores observados.
- Consideramos la función 
$$
L( \mathbf \beta) = \sum_{i=1}^n (y_i - \mu_i)^2 = 
   \sum_{i=1}^n  \bigg (y_i - \sum_{j=1}^p \beta_j x_{ij}\bigg )^2.
$$
- Los estimadores mínimo cuadráticos minimizan la función
$L( \mathbf \beta)$.
\item
$$
\hat{ \mathbf \beta} = ( \mathbf X^T  \mathbf X)^{-1}  \mathbf X^T  \mathbf y. 
$$

##  Predicciones y la matriz $H$  

- Las medias estimadas las obtenemos con
$$
\hat{ \mathbf \mu} =  \mathbf X ( \mathbf X^T  \mathbf X)^{-1}  \mathbf X^T  \mathbf y.
$$
- Si denotamos
$$
 \mathbf H =  \mathbf X ( \mathbf X^T  \mathbf X)^{-1}  \mathbf X^T,
$$
entonces
$$
\hat{ \mathbf \mu} =  \mathbf H  \mathbf y.
$$
\item
$$
E [\hat{ \mathbf \beta}] =  \mathbf \beta.
$$
\item
$$
var [\hat{ \mathbf \beta}] = \sigma^2 ( \mathbf X^T  \mathbf X)^{-1},
$$

	
##  Residuos y estimación de la varianza

- Los residuos son 
$$
\mathbf e =  \mathbf y -  \mathbf X \hat{ \mathbf \beta}.
$$ 
- Se tiene  
$$
E \sum_{i=1}^n \frac{ (Y_i - \hat{\mu}_i)^2}{n-p} = \sigma^2.
$$

# tamidata2::gse25171  

## Con `time`  y `Pi` 

- Analizamos los datos correspondientes a la sonda `261892_at`.

- Consideramos las variables fenotípicas `time` y `Pi`. 

```{r}
#| label: rls26
pacman::p_load(Biobase)
data(gse25171,package="tamidata2")
head(pData(gse25171),n=2)
```

```{r}
#| label: rls27
sel0 = which("261892_at"==fData(gse25171)[,"PROBEID"])
df0 = data.frame(pData(gse25171)[,c("time","Pi")],
                 expression=exprs(gse25171)[sel0,])
```

- Ajustamos un modelo con dos variables predictoras.  

```{r}
#| label: rls28
fit4 = lm(expression~ time + Pi,data=df0)
```
- La matriz modelo es 

```{r}
#| label: rls29
head(model.matrix(fit4),n=5)
```
- Los coeficientes los tenemos con 

```{r}
#| label: rls30
coef(fit4)
```

- Los residuos los tenemos con 

```{r}
#| label: rls31
resid(fit4)
```

-  El error estándar residual lo tenemos con 

```{r}
#| label: rls32
summary(fit4)$sigma
```  

## Con `time2`  y `Pi`   

- Analizamos los datos correspondientes a la sonda `261892_at`.

- Consideramos las variables fenotípicas `time2` y `Pi`. 

```{r}
#| label: rls33
pacman::p_load(Biobase)
data(gse25171,package="tamidata2")
head(pData(gse25171),n=2)
```

```{r}
#| label: rls34
sel0 = which("261892_at"==fData(gse25171)[,"PROBEID"])
df0 = data.frame(pData(gse25171)[,c("time2","Pi")],
                 expression=exprs(gse25171)[sel0,])
```
- Ajustamos un modelo con dos variables predictoras.

```{r}
#| label: rls35
fit5 = aov(expression~ time2 + Pi,data=df0)
```
- La matriz modelo es 

```{r}
#| label: rls36
head(model.matrix(fit5),n=5)
```
- La tabla de análisis de la varianza la tenemos con

```{r}
#| label: rls37
summary(aov(fit5))
```
- ¿Y un modelo con interacción?  

```{r}
#| label: rls38
fit6 = aov(expression~ time2 * Pi,data=df0)
```
- Vemos la tabla anova.  

```{r}
#| label: rls39
summary(fit6)
```
- No parece que tengamos interacción.  

# Muchos modelos 

## Un modelo por fila con `GSEAlm`
- Hemos elegido trabajar con una sonda elegida de un modo arbitrario.  
- Tenemos interés en todas las sondas.  
- Hemos de ajustar un modelo lineal para cada una de estas sondas.  
- Los predictores serán variables fenotípicas compartidas por todas las sondas.
- Los distintos modelos comparten los predictores y difieren en la variable respuesta.  

```{r}
#| label: rls40
pacman::p_load(GSEAlm)
data(gse25171,package="tamidata2")
fits1 = GSEAlm::lmPerGene(gse25171,~ time + Pi)
```
- La matriz de modelo común a todos los ajustes la tenemos con
```{r}
#| label: rls41
head(fits1$x)
```
- Los coeficientes de cada ajuste con 
```{r}
#| label: rls42
#| results: false  
fits1$coefficients
```
-  Para la primera sonda tenemos los coeficientes
```{r}
#| label: rls42b
fits1$coefficients[,1]
```
- Las varianzas del error estimadas son 
```{r}
#| label: rls43
#| results: false
fits1$sigmaSqr
```
Y los estadísticos correspondientes a los tests de coeficientes nulos con
```{r}
#| label: rls44
#| results: false
fits1$tstat
``` 

## Un modelo por fila con `limma`  

```{r}
#| label: rls45
pacman::p_load(limma)
```
- Ajustamos los modelos. 

```{r}
#| label: rls46
design = model.matrix(~ pData(gse25171)[,"time"] +  pData(gse25171)[,"Pi"])
fits2 = limma::lmFit(gse25171,design)
```  
- Los coeficientes los tenemos con 

```{r}
#| label: rls47
head(fits2$coefficients,n=2)
```

- Los coeficientes del ajuste para la primera sonda son
```{r}
#| label: rls48
fits2$coefficients[1,] 
```

 
