Análisis de datos PISA, los valores plausibles

En este post puedes descargar los ejemplos de código R para trabajar con los valores plausibles en la base de datos de PISA, para calcular medias, diferencias de medias o regresiones lineales de las calificaciones de los alumnos, usando los pesos replicados para obtener los errores estándar.

Este post está relacionado con el artículo cálculos con valores plausibles en la base de datos PISA.

En este enlace podéis descargar el programa R para Windows.

En este enlace podéis descargar el código R de ejemplo para cálculo con valores plausibles. Estas funciones trabajan con un dataframe del que previamente habremos eliminado las filas con valores faltantes, por simplicidad. En lo que sigue vamos a hacer un ligero repaso de cada una de estas funciones y sus parámetros.

Ejemplos de cálculo de la media y la desviación típica con valores plausibles

En el script tenemos dos funciones para calcular la media y la desviación típica de los valores plausibles en un conjunto de datos, junto con sus errores estándar, calculados a partir de los pesos replicados, como vimos en el artículo cálculo de errores estándar con pesos replicados en la base de datos PISA.

Ejemplo 1, cálculo de la media y la desviación típica

En este ejemplo, calculamos el valor correspondiente a la media y la desviación típica, junto con sus errores típicos, para un conjunto de valores plausibles.

La función es wght_meansd_pv, y este es el código:

wght_meansd_pv<-function(sdata,pv,wght,brr) {
mmeans<-c(0, 0, 0, 0);
mmeanspv<-rep(0,length(pv));
stdspv<-rep(0,length(pv));
mmeansbr<-rep(0,length(pv));
stdsbr<-rep(0,length(pv));
names(mmeans)<-c("MEAN","SE-MEAN","STDEV","SE-STDEV");
swght<-sum(sdata[,wght]);
for (i in 1:length(pv)) {
mmeanspv[i]<-sum(sdata[,wght]*sdata[,pv[i]])/swght;
stdspv[i]<-sqrt((sum(sdata[,wght]*(sdata[,pv[i]]^2))/swght)-
mmeanspv[i]^2);
for (j in 1:length(brr)) {
sbrr<-sum(sdata[,brr[j]]);
mbrrj<-sum(sdata[,brr[j]]*sdata[,pv[i]])/sbrr;
mmeansbr[i]<-mmeansbr[i] + (mbrrj - mmeanspv[i])^2;
stdsbr[i]<-stdsbr[i] +
(sqrt((sum(sdata[,brr[j]]*(sdata[,pv[i]]^2))/sbrr)-mbrrj^2) -
stdspv[i])^2;
}
}
mmeans[1]<-sum(mmeanspv) / length(pv);
mmeans[2]<-sum((mmeansbr * 4) / length(brr)) / length(pv);
mmeans[3]<-sum(stdspv) / length(pv);
mmeans[4]<-sum((stdsbr * 4) / length(brr)) / length(pv);
ivar <- c(0,0);
for (i in 1:length(pv)) {
ivar[1] <- ivar[1] + (mmeanspv[i] - mmeans[1])^2;
ivar[2] <- ivar[2] + (stdspv[i] - mmeans[3])^2;
}
ivar = (1 + (1 / length(pv))) * (ivar / (length(pv) - 1));
mmeans[2]<-sqrt(mmeans[2] + ivar[1]);
mmeans[4]<-sqrt(mmeans[4] + ivar[2]);
return(mmeans);
}

Para calcular la media y la desviación típica, tenemos que sumar cada uno de los cinco valores plausibles multiplicado por el peso del alumno, y, a continuación, realizar la media de los resultados parciales de cada valor. Para calcular el error típico usamos el método de los pesos replicados, pero tenemos que añadir la varianza de imputación, entre los cinco valores plausibles, lo cual hacemos mediante la variable ivar.

En el parámetro sdata pasaremos el dataframe con los datos. Los nombres o índices de columna de los valores plausibles se pasan en un vector en el parámetro pv, mientras que los parámetros wght (índice o nombre de la columna con el peso del alumno) y brr (vector de índices o nombres de columna para los pesos replicados) se usan como hemos visto en anteriores artículos.

Como resultado obtenemos un vector con cuatro posiciones, la primera para la nota media, la segunda para su error estándar, la tercera para la desviación típica y la cuarta para el error estándar de la desviación típica.

Ejemplo 2, cálculo de la media y la desviación típica agrupando por los niveles de una serie de factores

En este ejemplo se realiza el mismo cálculo del ejemplo anterior, pero esta vez agrupando por los niveles de una o más columnas con datos de tipo factor, como puede ser el sexo del alumno o el grado en el que estaba en el momento del examen.

La función es wght_meansdfact_pv, y el código es el siguiente:

wght_meansdfact_pv<-function(sdata,pv,cfact,wght,brr) {
nc<-0;
for (i in 1:length(cfact)) {
nc <- nc + length(levels(as.factor(sdata[,cfact[i]])));
}
mmeans<-matrix(ncol=nc,nrow=4);
mmeans[,]<-0;
cn<-c();
for (i in 1:length(cfact)) {
for (j in 1:length(levels(as.factor(sdata[,cfact[i]])))) {
cn<-c(cn, paste(names(sdata)[cfact[i]],
levels(as.factor(sdata[,cfact[i]]))[j],sep="-"));
}
}
colnames(mmeans)<-cn;
rownames(mmeans)<-c("MEAN","SE-MEAN","STDEV","SE-STDEV");
ic<-1;
for(f in 1:length(cfact)) {
for (l in 1:length(levels(as.factor(sdata[,cfact[f]])))) {
rfact<-sdata[,cfact[f]]==levels(as.factor(sdata[,cfact[f]]))[l];
swght<-sum(sdata[rfact,wght]);
mmeanspv<-rep(0,length(pv));
stdspv<-rep(0,length(pv));
mmeansbr<-rep(0,length(pv));
stdsbr<-rep(0,length(pv));
for (i in 1:length(pv)) {
mmeanspv[i]<-sum(sdata[rfact,wght]*sdata[rfact,pv[i]])/swght;
stdspv[i]<-sqrt((sum(sdata[rfact,wght] *
(sdata[rfact,pv[i]]^2))/swght)-mmeanspv[i]^2);
for (j in 1:length(brr)) {
sbrr<-sum(sdata[rfact,brr[j]]);
mbrrj<-sum(sdata[rfact,brr[j]]*sdata[rfact,pv[i]])/sbrr;
mmeansbr[i]<-mmeansbr[i] + (mbrrj - mmeanspv[i])^2;
stdsbr[i]<-stdsbr[i] + (sqrt((sum(sdata[rfact,brr[j]] *
(sdata[rfact,pv[i]]^2))/sbrr)-mbrrj^2) - stdspv[i])^2;
}
}
mmeans[1, ic]<- sum(mmeanspv) / length(pv);
mmeans[2, ic]<-sum((mmeansbr * 4) / length(brr)) / length(pv);
mmeans[3, ic]<- sum(stdspv) / length(pv);
mmeans[4, ic]<-sum((stdsbr * 4) / length(brr)) / length(pv);
ivar <- c(sum((mmeanspv - mmeans[1, ic])^2), sum((stdspv -
mmeans[3, ic])^2));
ivar = (1 + (1 / length(pv))) * (ivar / (length(pv) - 1));
mmeans[2, ic]<-sqrt(mmeans[2, ic] + ivar[1]);
mmeans[4, ic]<-sqrt(mmeans[4, ic] + ivar[2]);
ic<-ic + 1;
}
}
return(mmeans);
}

Además de los parámetros que ya hemos visto en el ejemplo anterior, con el mismo uso y significado, tenemos el parámetro cfact, en el que deberemos pasar un vector con los índices o nombres de columna de los factores por cuyos niveles queremos agrupar los datos.

El resultado se devuelve en una matriz con cuatro filas, la primera para las medias, la segunda para sus errores estándar, la tercera para la desviación típica y la cuarta para el error estándar de la desviación típica. En cada columna tendremos el valor correspondiente para cada uno de los niveles de cada uno de los factores.

Ejemplos de cálculo de diferencias de medias

En los dos ejemplos que siguen, veremos cómo calcular diferencias de medias de valores plausibles y sus errores típicos usando los pesos replicados.

Ejemplo 3, cálculo de la diferencia de medias entre países

Esta función trabaja sobre un dataframe que contiene datos para varios países, y calcula la diferencia de medias entre cada pareja de dos países.

La función es wght_meandiffcnt_pv, y el código es el siguiente:

wght_meandiffcnt_pv<-function(sdata,pv,cnt,wght,brr) {
nc<-0;
for (j in 1:(length(levels(as.factor(sdata[,cnt])))-1)) {
for(k in (j+1):length(levels(as.factor(sdata[,cnt])))) {
nc <- nc + 1;
}
}
mmeans<-matrix(ncol=nc,nrow=2);
mmeans[,]<-0;
cn<-c();
for (j in 1:(length(levels(as.factor(sdata[,cnt])))-1)) {
for(k in (j+1):length(levels(as.factor(sdata[,cnt])))) {
cn<-c(cn, paste(levels(as.factor(sdata[,cnt]))[j],
levels(as.factor(sdata[,cnt]))[k],sep="-"));
}
}
colnames(mmeans)<-cn;
rn<-c("MEANDIFF", "SE");
rownames(mmeans)<-rn;
ic<-1;
for (l in 1:(length(levels(as.factor(sdata[,cnt])))-1)) {
for(k in (l+1):length(levels(as.factor(sdata[,cnt])))) {
rcnt1<-sdata[,cnt]==levels(as.factor(sdata[,cnt]))[l];
rcnt2<-sdata[,cnt]==levels(as.factor(sdata[,cnt]))[k];
swght1<-sum(sdata[rcnt1,wght]);
swght2<-sum(sdata[rcnt2,wght]);
mmeanspv<-rep(0,length(pv));
mmcnt1<-rep(0,length(pv));
mmcnt2<-rep(0,length(pv));
mmeansbr1<-rep(0,length(pv));
mmeansbr2<-rep(0,length(pv));
for (i in 1:length(pv)) {
mmcnt1<-sum(sdata[rcnt1,wght]*sdata[rcnt1,pv[i]])/swght1;
mmcnt2<-sum(sdata[rcnt2,wght]*sdata[rcnt2,pv[i]])/swght2;
mmeanspv[i]<- mmcnt1 - mmcnt2;
for (j in 1:length(brr)) {
sbrr1<-sum(sdata[rcnt1,brr[j]]);
sbrr2<-sum(sdata[rcnt2,brr[j]]);
mmbrj1<-sum(sdata[rcnt1,brr[j]]*sdata[rcnt1,pv[i]])/sbrr1;
mmbrj2<-sum(sdata[rcnt2,brr[j]]*sdata[rcnt2,pv[i]])/sbrr2;
mmeansbr1[i]<-mmeansbr1[i] + (mmbrj1 - mmcnt1)^2;
mmeansbr2[i]<-mmeansbr2[i] + (mmbrj2 - mmcnt2)^2;
}
}
mmeans[1,ic]<-sum(mmeanspv) / length(pv);
mmeansbr1<-sum((mmeansbr1 * 4) / length(brr)) / length(pv);
mmeansbr2<-sum((mmeansbr2 * 4) / length(brr)) / length(pv);
mmeans[2,ic]<-sqrt(mmeansbr1^2 + mmeansbr2^2);
ivar <- 0;
for (i in 1:length(pv)) {
ivar <- ivar + (mmeanspv[i] - mmeans[1,ic])^2;
}
ivar = (1 + (1 / length(pv))) * (ivar / (length(pv) - 1));
mmeans[2,ic]<-sqrt(mmeans[2,ic] + ivar);
ic<-ic + 1;
}
}
return(mmeans);
}

Aquí el cálculo de los errores estándar es diferente. Se debe calcular el error típico para cada país por separado, y después obtener la raíz cuadrada de la suma de los dos cuadrados, pues los datos de cada país son independientes del resto.

De nuevo los parámetros son los mismos que en las funciones anteriores. Tenemos el nuevo parámetro cnt, en el cual deberemos pasar el índice o el nombre de la columna con el país.

El resultado es una matriz con dos filas, la primera para las diferencias y la segunda para los errores estándar, y una columna para la diferencia entre cada una de las combinaciones de los países.

Ejemplo 4, cálculo de las diferencias de medias agrupando por los niveles de una serie de factores

Con esta función agrupamos por los niveles de una serie de factores y calculamos también las diferencias de medias dentro de cada país, además de las diferencias de medias entre países.

La función es wght_meandifffactcnt_pv, y el código es el siguiente:

wght_meandifffactcnt_pv<-function(sdata,pv,cnt,cfact,wght,brr) {
lcntrs<-vector('list',1 + length(levels(as.factor(sdata[,cnt]))));
for (p in 1:length(levels(as.factor(sdata[,cnt])))) {
names(lcntrs)[p]<-levels(as.factor(sdata[,cnt]))[p];
}
names(lcntrs)[1 + length(levels(as.factor(sdata[,cnt])))]<-"BTWNCNT";
nc<-0;
for (i in 1:length(cfact)) {
for (j in 1:(length(levels(as.factor(sdata[,cfact[i]])))-1)) {
for(k in (j+1):length(levels(as.factor(sdata[,cfact[i]])))) {
nc <- nc + 1;
}
}
}
cn<-c();
for (i in 1:length(cfact)) {
for (j in 1:(length(levels(as.factor(sdata[,cfact[i]])))-1)) {
for(k in (j+1):length(levels(as.factor(sdata[,cfact[i]])))) {
cn<-c(cn, paste(names(sdata)[cfact[i]],
levels(as.factor(sdata[,cfact[i]]))[j],
levels(as.factor(sdata[,cfact[i]]))[k],sep="-"));
}
}
}
rn<-c("MEANDIFF", "SE");
for (p in 1:length(levels(as.factor(sdata[,cnt])))) {
mmeans<-matrix(ncol=nc,nrow=2);
mmeans[,]<-0;
colnames(mmeans)<-cn;
rownames(mmeans)<-rn;
ic<-1;
for(f in 1:length(cfact)) {
for (l in 1:(length(levels(as.factor(sdata[,cfact[f]])))-1)) {
for(k in (l+1):length(levels(as.factor(sdata[,cfact[f]])))) {
rfact1<-
(sdata[,cfact[f]] ==
levels(as.factor(sdata[,cfact[f]]))[l]) &
(sdata[,cnt]==levels(as.factor(sdata[,cnt]))[p]);
rfact2<-
(sdata[,cfact[f]] ==
levels(as.factor(sdata[,cfact[f]]))[k]) &
(sdata[,cnt]==levels(as.factor(sdata[,cnt]))[p]);
swght1<-sum(sdata[rfact1,wght]);
swght2<-sum(sdata[rfact2,wght]);
mmeanspv<-rep(0,length(pv));
mmeansbr<-rep(0,length(pv));
for (i in 1:length(pv)) {
mmeanspv[i]<-(sum(sdata[rfact1,wght] *
sdata[rfact1,pv[i]])/swght1) -
(sum(sdata[rfact2,wght] *
sdata[rfact2,pv[i]])/swght2);
for (j in 1:length(brr)) {
sbrr1<-sum(sdata[rfact1,brr[j]]);
sbrr2<-sum(sdata[rfact2,brr[j]]);
mmbrj<-(sum(sdata[rfact1,brr[j]] *
sdata[rfact1,pv[i]])/sbrr1) -
(sum(sdata[rfact2,brr[j]] *
sdata[rfact2,pv[i]])/sbrr2);
mmeansbr[i]<-mmeansbr[i] + (mmbrj - mmeanspv[i])^2;
}
}
mmeans[1,ic]<-sum(mmeanspv) / length(pv);
mmeans[2,ic]<-sum((mmeansbr * 4) / length(brr)) /
length(pv);
ivar <- 0;
for (i in 1:length(pv)) {
ivar <- ivar + (mmeanspv[i] - mmeans[1,ic])^2;
}
ivar = (1 + (1 / length(pv))) * (ivar / (length(pv) - 1));
mmeans[2,ic]<-sqrt(mmeans[2,ic] + ivar);
ic<-ic + 1;
}
}
}
lcntrs[[p]]<-mmeans;
}
pn<-c();
for (p in 1:(length(levels(as.factor(sdata[,cnt])))-1)) {
for (p2 in (p + 1):length(levels(as.factor(sdata[,cnt])))) {
pn<-c(pn,
paste(levels(as.factor(sdata[,cnt]))[p],
levels(as.factor(sdata[,cnt]))[p2],sep="-"));
}
}
mbtwmeans<-array(0, c(length(rn), length(cn), length(pn)));
nm <- vector('list',3);
nm[[1]]<-rn;
nm[[2]]<-cn;
nm[[3]]<-pn;
dimnames(mbtwmeans)<-nm;
pc<-1;
for (p in 1:(length(levels(as.factor(sdata[,cnt])))-1)) {
for (p2 in (p + 1):length(levels(as.factor(sdata[,cnt])))) {
ic<-1;
for(f in 1:length(cfact)) {
for (l in 1:(length(levels(as.factor(sdata[,cfact[f]])))-1)) {
for(k in (l+1):length(levels(as.factor(sdata[,cfact[f]]))))
{
mbtwmeans[1,ic,pc]<-lcntrs[[p]][1,ic] -
lcntrs[[p2]][1,ic];
mbtwmeans[2,ic,pc]<-sqrt((lcntrs[[p]][2,ic]^2) +
(lcntrs[[p2]][2,ic]^2));
ic<-ic + 1;
}
}
}
pc<-pc+1;
}
}
lcntrs[[1 + length(levels(as.factor(sdata[,cnt])))]]<-mbtwmeans;
return(lcntrs);
}

A los parámetros de la función del ejemplo anterior, hemos añadido cfact, donde se pasa el vector con los índices o nombres de columna de los factores.

En este caso, los datos se devuelven en una lista. Para cada uno de los países existe un elemento de la lista que contienen una matriz con dos filas, una para las diferencias y otra para los errores estándar, y una columna para cada combinación posible de dos niveles de cada uno de los factores, a partir de la cual se calculan las diferencias.

En el último elemento de la lista, se devuelve un array de tres dimensiones, una de ellas para cada combinación de dos países, y las otras dos forman una matriz con la misma estructura de filas y columnas que las matrices de las posiciones correspondientes a cada país.

Regresiones lineales con valores plausibles

En este último ejemplo, veremos un función para realizar regresiones lineales en las que las variables dependientes son los valores plausibles, obteniendo los coeficientes de regresión y sus errores estándar.

Ejemplo 5, regresión lineal usando los valores plausibles como variables dependientes

La función es wght_lmpv, y este es el código:

wght_lmpv<-function(sdata,frml,pv,wght,brr) {
listlm <- vector('list', 2 + length(pv));
listbr <- vector('list', length(pv));
for (i in 1:length(pv)) {
if (is.numeric(pv[i])) {
names(listlm)[i] <- colnames(sdata)[pv[i]];
frmlpv <- as.formula(paste(colnames(sdata)[pv[i]],frml,sep="~"));
}
else {
names(listlm)[i]<-pv[i];
frmlpv <- as.formula(paste(pv[i],frml,sep="~"));
}
listlm[[i]] <- lm(frmlpv, data=sdata, weights=sdata[,wght]);
listbr[[i]] <- rep(0,2 + length(listlm[[i]]$coefficients));
for (j in 1:length(brr)) {
lmb <- lm(frmlpv, data=sdata, weights=sdata[,brr[j]]);
listbr[[i]]<-listbr[[i]] + c((listlm[[i]]$coefficients -
lmb$coefficients)^2,(summary(listlm[[i]])$r.squared-
summary(lmb)$r.squared)^2,(summary(listlm[[i]])$adj.r.squared-
summary(lmb)$adj.r.squared)^2);
}
listbr[[i]] <- (listbr[[i]] * 4) / length(brr);
}
cf <- c(listlm[[1]]$coefficients,0,0);
names(cf)[length(cf)-1]<-"R2";
names(cf)[length(cf)]<-"ADJ.R2";
for (i in 1:length(cf)) {
cf[i] <- 0;
}
for (i in 1:length(pv)) {
cf<-(cf + c(listlm[[i]]$coefficients, summary(listlm[[i]])$r.squared,
summary(listlm[[i]])$adj.r.squared));
}
names(listlm)[1 + length(pv)]<-"RESULT";
listlm[[1 + length(pv)]]<- cf / length(pv);
names(listlm)[2 + length(pv)]<-"SE";
listlm[[2 + length(pv)]] <- rep(0, length(cf));
names(listlm[[2 + length(pv)]])<-names(cf);
for (i in 1:length(pv)) {
listlm[[2 + length(pv)]] <- listlm[[2 + length(pv)]] + listbr[[i]];
}
ivar <- rep(0,length(cf));
for (i in 1:length(pv)) {
ivar <- ivar + c((listlm[[i]]$coefficients - listlm[[1 +
length(pv)]][1:(length(cf)-2)])^2,(summary(listlm[[i]])$r.squared -
listlm[[1 + length(pv)]][length(cf)-1])^2,
(summary(listlm[[i]])$adj.r.squared - listlm[[1 +
length(pv)]][length(cf)])^2);
}
ivar = (1 + (1 / length(pv))) * (ivar / (length(pv) - 1));
listlm[[2 + length(pv)]] <- sqrt((listlm[[2 + length(pv)]] / length(pv)) +
ivar);
return(listlm);
}

En esta función, debemos pasar la parte derecha de la fórmula como una cadena de texto en el parámetro frml, por ejemplo, si las variables independientes fuesen HISEI y ST03Q01, pasaríamos la cadena de texto “HISEI+ST03Q01”. La función calculará un modelo lineal para cada uno de los valores plausibles con la función lm y, a partir de estos, construirá el modelo final y calculará los errores estándar.

Como resultado obtendremos una lista, con una posición para los coeficientes de cada uno de los modelos de cada valor plausible, otra con los coeficientes del resultado final, y otra más con los errores estándar correspondientes a dichos coeficientes.