[R-es] diseño matriz problema model.matrix

Javier Marcuzzi j@v|er@ruben@m@rcuzz| @end|ng |rom gm@||@com
Lun Ene 27 04:06:04 CET 2025


Estimados

Copio y pego un código reproducible donde está la pregunta, yo lo realicé con en un archivo cmd pero si no me equivoco en rstudio o cualquier otro debería ser leído sin problemas al cortar y pegar.


## Consulta y descripción del problema

De un libro estoy siguiendo unos problemas, este libro se está convirtiendo en referencia en un área muy específica, tal es así que el año pasado realizaron la cuarta edición del mismo.

Como se va convirtiendo en referencia y ser un libro donde no utilizan programas o lenguajes informáticos, distintas personas realizan en diversos lenguajes replicas de los ejemplos del libro. Tal es así, que la misma editorial publica en R la respuestas aportada por terceros.

En R tengo tres versiónes, la mía que da una falla, la de un tercero citada por la editorial y la de un Phd de una universidad de Estados Unidos.

La diferencia entre el Phd y el recomendado por la editorial, es que el primero utiliza una función para diseñar la matriz, función casi idéntica a otra universidad de Canadá, en cambio el recomendado utiliza model.matrix. Ambas formas son idénticas pero luego difieren, y esta diferencia la colocaré en ejemplos en esta consulta a la lista.

El libro tiene todo escrito en matemática, álgebra, me refiero a que no tiene ni una línea de código, por lo que cada persona es libre de utilizar papel y lápiz o la computadora.

El ejemplo en R con la función escrita por el phd de Estados Unidos, coincide con el libro, es fiel a este. Sin embargo continuando con los pasos para la resolución del problema, hay error, el mensaje es el siguiente: Error in solve.default(LHS, RHS) : Lapack routine dgesv: system is exactly singular: U\[21,21\] = 0

En cambio el otro ejemplo recomendado por la editorial, no da este error y llega a la resolución, pero no coincide con el libro, no es una copia fiel en el "armado" de la matriz.

La mejor forma de explicar la pregunta es compartiendo un código reproducible.

## Datos

Los datos son los mismos, los copio antes de continuar el trabajo, uno se llama data y el otro data0. La copia es porque al colocar "factor" no puedo utilizarlos en ambos casos. Se comprenderá en el ejemplo.

```{r}
cow = rep(seq(4, 8), each = 2)
fy = c(201, 280, 150, 200, 160, 190, 180, 250, 285, 300)
hys = c(1, 3, 1, 4, 2, 3, 1, 3, 2, 4)
pairity = rep(c(1,2), 5)

data = data.frame(cow, pairity, hys, fy)
data0 = data.frame(cow, pairity, hys, fy)
```

## Preparación de datos

Para utilizar con model.matrix, utilizo factor. En este caso se llaman data.

```{r}
data$pairity = as.factor(pairity)
data$hys = as.factor(hys)
data$hys = relevel(factor(hys), ref = "1")
data$cow = factor(data$cow)
```

Para utilizar la función de diseño, uso los datos originales (sin factor). Como no tienen "trabajo extra" para preprar los datos se llaman datos0.

## La función de diseño de matriz es la siguiente:

```{r}
# función para crear el diseño de matriz 
desgn <- function(v) {
  if (is.numeric(v)) {
    va = v
    mrow = length(va)
    mcol = max(va)
  }
  if (is.character(v)) {
    vf = factor(v)
    # Guarda el indice para cada nivel
    va = as.numeric(vf)
    mrow = length(va)
    mcol = length(levels(vf))
  }
  
  # Inicio la martiz
  X = matrix(data = c(0), nrow = mrow, ncol = mcol)
  
  for (i in 1:mrow) {
    ic = va[i]
    X[i, ic] = 1
  }
  return(X)
}
```

## Armado de las matrices

### Forma 1

Utilizando la función de diseño, hace una fiel coincidencia con el libro, me refiero a que el libro tiene todo escrito en forma matemática, sin utilizar lenguajes informáticos, se puede hacer con lápiz y papel, no hay otro requerimiento. Pero luego este armado de matrices da error (singularidad).

Como se puede reproducir, el ejemplo son dos matrices que se unen con cbind.

```{r}
# Matriz X, relaciona los registros de HYS y partos
# dideño para HYS
x1 = desgn(data0$hys)
x1
# diseño para Partos
x2 = desgn(data0$pairity)
x2
# Matriz X segín libro
XLibro = cbind(x1,x2)

t(XLibro)
```


Más información sobre la lista de distribución R-help-es