Análisis de datos PISA, Cálculo de los errores estándar con los pesos replicados

En este post puedes descargar los ejemplos de código R sobre el cálculo de los errores estándar de la media, desviación típica, proporciones o diferencias de medias en los datos de la base de datos de PISA, utilizando el método de los pesos replicados.

Este post está relacionado con el artículo cálculo de errores estándar con pesos replicados en la base de datos PISA.

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

Todos los ejemplos se deben realizar sobre un dataframe del que se hayan eliminado los valores faltantes.

Ejemplos de cálculo del error típico de la media y la desviación típica

En este enlace podéis descargar los ejemplos de código para calcular la media y la desviación típica. Se trata de tres funciones que operan sobre un dataframe. Con una de ellas calculamos la media y la desviación típica de una serie de columnas, con otra realizamos los mismos cálculos, pero sobre un dataframe que contiene más de un país y con la tercera agrupamos los datos por los diferentes niveles de una serie de factores antes de realizar los cálculos.

Ejemplo 1, cálculo de la media, la desviación típica y sus correspondientes errores estándar

La primera función se llama wght_meansd, y este es el código:

wght_meansd<-function(sdata,cmean,wght,brr) {
mmeans<-matrix(ncol=length(cmean),nrow=4);
mmeans[,]<-0;
if (is.numeric(cmean)) {
colnames(mmeans)<-names(sdata)[cmean];
}
else {
colnames(mmeans)<-cmean;
}
rownames(mmeans)<-c("MEAN","SE-MEAN","STDEV","SE-STDEV");
swght<-sum(sdata[,wght]);
for (i in 1:length(cmean)) {
mmeans[1,i]<-sum(sdata[,wght]*sdata[,cmean[i]])/swght;
mmeans[3,i]<-sqrt((sum(sdata[,wght]*(sdata[,cmean[i]]^2))/swght)-
mmeans[1,i]^2);
for (j in 1:length(brr)) {
sbrr<-sum(sdata[,brr[j]]);
mmbrj<-sum(sdata[,brr[j]]*sdata[,cmean[i]])/sbrr;
mmeans[2,i]<-mmeans[2,i] + (mmbrj - mmeans[1,i])^2;
mmeans[4,i]<-mmeans[4,i] +
(sqrt((sum(sdata[,brr[j]]*(sdata[,cmean[i]]^2))/sbrr)-mmbrj^2) -
mmeans[3,i])^2;
}
}
mmeans[2,]<-sqrt((mmeans[2,] * 4) / length(brr));
mmeans[4,]<-sqrt((mmeans[4,] * 4) / length(brr));
return(mmeans);
}

En el parámetro sdata pasaremos el dataframe con los datos. cmean es un vector con los índices o los nombres de las columnas de las que queremos calcular la media y la desviación típica. wght es el índice o el nombre de la columna con el peso del alumno (W_FSTUWT), y brr es un vector con los índices o los nombres de las columnas con los 80 pesos replicados (W_FSTR1 a W_FSTR80).

Para el cálculo de la media y la desviación típica, los valores en cada columna se deben ponderar usando los pesos de los alumnos. Los errores típicos se obtienen calculando la media o la desviación típica usando cada uno de los pesos replicados y sumando el cuadrado de sus diferencias con la media calculada usando el peso del alumno. El resultado final de todas estas sumas se divide por 20 y se calcula la raíz cuadrada.

El resultado se devuelve en una matriz con las mismas columnas que las pasadas en el parámetro cmean, y con cuatro filas, las dos primeras corresponden a las medias y su error típico, y las dos últimas a las desviaciones típicas y sus errores estándar.

Ejemplo 2, el mismo cálculo anterior, pero aplicado a una muestra con varios países

La función de este ejemplo es wght_meansdbycnt, y este es el código:

wght_meansdbycnt<-function(sdata,cmean,cnt,wght,brr) {
lcntrs<-vector('list',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];
}
for (p in 1:length(levels(as.factor(sdata[,cnt])))) {
mmeans<-matrix(ncol=length(cmean),nrow=4);
mmeans[,]<-0;
if (is.numeric(cmean)) {
colnames(mmeans)<-names(sdata)[cmean];
}
else {
colnames(mmeans)<-cmean;
}
rownames(mmeans)<-c("MEAN","SE-MEAN","STDEV","SE-STDEV");
rcnt <- sdata[,cnt]==levels(as.factor(sdata[,cnt]))[p];
swght<-sum(sdata[rcnt,wght]);
for (i in 1:length(cmean)) {
mmeans[1,i]<-sum(sdata[rcnt,wght]*sdata[rcnt,cmean[i]])/swght;
mmeans[3,i]<-sqrt((sum(sdata[,wght]*(sdata[,cmean[i]]^2))/swght)-
mmeans[1,i]^2);
for (j in 1:length(brr)) {
sbrr<-sum(sdata[rcnt,brr[j]]);
mmbrj<-sum(sdata[rcnt,brr[j]]*sdata[rcnt,cmean[i]])/sbrr;
mmeans[2,i]<-mmeans[2,i] + (mmbrj - mmeans[1,i])^2;
mmeans[4,i]<-mmeans[4,i] +
(sqrt((sum(sdata[rcnt,brr[j]]*(sdata[rcnt,cmean[i]]^2))/sbrr)-
mmbrj^2) - mmeans[3,i])^2;
}
}
mmeans[2,]<-sqrt((mmeans[2,] * 4) / length(brr));
mmeans[4,]<-sqrt((mmeans[4,] * 4) / length(brr));
lcntrs[[p]]<-mmeans;
}
return(lcntrs);
}

La función tiene un parámetro adicional, cnt, en el que se pasa el índice o el nombre de la columna con el país. Los cálculos realizados son los mismos que en el caso anterior, pero el resultado se devuelve en una lista de matrices, cada una correspondiente a uno de los países.

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

En este ejemplo realizamos el cálculo de la media y desviación típica, junto con sus errores estándar, pero separándolos según los niveles de una serie de factores, como pueden ser el sexo o el nivel de estudios de los padres.

La función es wght_meansdfacty este es el código:

wght_meansdfact<-function(sdata,cmean,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 * length(cmean));
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;
rn<-c();
for (i in 1:length(cmean)) {
rn<-c(rn, paste(names(sdata)[cmean[i]],"MEAN",sep="-"));
rn<-c(rn, paste(names(sdata)[cmean[i]],"MEAN","SE",sep="-"));
rn<-c(rn, paste(names(sdata)[cmean[i]],"STDEV",sep="-"));
rn<-c(rn, paste(names(sdata)[cmean[i]],"STDEV","SE",sep="-"));
}
rownames(mmeans)<-rn;
ic<-1;
for(f in 1:length(cfact)) {
for (l in 1:length(levels(as.factor(sdata[,cfact[f]])))) {
ir<-1;
rfact<-sdata[,cfact[f]]==levels(as.factor(sdata[,cfact[f]]))[l];
swght<-sum(sdata[rfact,wght]);
for (i in 1:length(cmean)) {
mmeans[ir,ic]<-sum(sdata[rfact,wght] *
sdata[rfact,cmean[i]])/swght;
mmeans[ir+2,ic]<-sqrt((sum(sdata[rfact,wght] *
(sdata[rfact,cmean[i]]^2))/swght)-mmeans[ir,ic]^2);
for (j in 1:length(brr)) {
sbrr<-sum(sdata[rfact,brr[j]]);
mmbrj<-sum(sdata[rfact,brr[j]]*sdata[rfact,cmean[i]])/sbrr;
mmeans[ir+1,ic]<-mmeans[ir+1,ic] + (mmbrj -
mmeans[ir,ic])^2;
mmeans[ir+3,ic]<-mmeans[ir+3,ic] +
(sqrt((sum(sdata[rfact,brr[j]] *
(sdata[rfact,cmean[i]]^2))/sbrr)-mmbrj^2) -
mmeans[ir+2,ic])^2;
}
mmeans[ir+1,ic]<-sqrt((mmeans[ir+1,ic] * 4) / length(brr));
mmeans[ir+3,ic]<-sqrt((mmeans[ir+3,ic] * 4) / length(brr));
ir<-ir+4;
}
ic<-ic + 1;
}
}
return(mmeans);
}

Esta vez, la agrupación la realizamos según los diferentes niveles de los factores cuyo índice o nombre de columna se pasa en un vector en el parámetro cfact. Como resultado obtenemos una matriz similar a la del ejemplo 1, pero en cada columna tendremos la media y desviación típica para un determinado nivel de factor. Por ejemplo, tendremos la media de alumnos varones en una columna y en otra las mujeres.

Ejemplos de cálculo de proporciones y su error típico

Las funciones de estos ejemplos realizan un recuento de niveles de factor para calcular sus proporciones en una muestra, así como el error típico de dichas proporciones. En este enlace puedes descargar los ejemplos de código para el cálculo de proporciones.

En estos y sucesivos ejemplos, los parámetros wght y brr tienen el mismo uso y significado que en el ejemplo 1 del punto anterior.

Ejemplo 4, cálculo simple de las proporciones de una serie de factores en una muestra

En este ejemplo simplemente obtendremos el porcentaje de cada uno de los niveles de una serie de factores de los casos que se dan en la muestra.

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

wght_ppc<-function(sdata,cppc,wght,brr) {
nc<-0;
for (i in 1:length(cppc)) {
nc <- nc + length(levels(as.factor(sdata[,cppc[i]])));
}
mppc<-matrix(ncol=nc,nrow=2);
mppc[,]<-0;
cn<-c();
for (i in 1:length(cppc)) {
for (j in 1:length(levels(as.factor(sdata[,cppc[i]])))) {
cn<-c(cn, paste(names(sdata)[cppc[i]],
levels(as.factor(sdata[,cppc[i]]))[j], sep="-"));
}
}
colnames(mppc)<-cn;
rownames(mppc)<-c("PPC","SE");
swght<-sum(sdata[,wght]);
ix<-1;
for (i in 1:length(cppc)) {
for (j in 1:length(levels(as.factor(sdata[,cppc[i]])))) {
rfact<-sdata[,cppc[i]]==levels(as.factor(sdata[,cppc[i]]))[j];
mppc[1,ix]<-sum(sdata[rfact,wght]) / swght;
for (k in 1:length(brr)) {
sbrr<-sum(sdata[,brr[k]]);
ppcbrr<-sum(sdata[rfact,brr[k]]) / sbrr;
mppc[2,ix]<-mppc[2,ix] + (ppcbrr-mppc[1,ix])^2;
}
x<-ix + 1;
}
}
mppc[2,]<-sqrt((mppc[2,] * 4) / length(brr));
return(mppc);
}

Para realizar el recuento, en realidad lo que hacemos es sumar los pesos en lugar de contar los casos. Para el cálculo de los errores estándar empelamos un procedimiento similar al del caso de las medias: calculamos la proporción para cada uno de los pesos replicados y calculamos el error estándar por el cuadrado de las diferencias con la proporción obtenida con el peso del alumno.

El parámetro sdata es para el dataframe con los datos, y cppc es un vector con los índices o nombres de las columnas con los factores.

Como resultado obtendremos una matriz con dos filas. En la primera tendremos las proporciones, entre 0 y 1, de los diferentes niveles de los factores, una columna por nivel. En la segunda fila estarán los errores estándar correspondientes.

Ejemplo 5, cálculo de las proporciones de una serie de factores, combinando sus distintos niveles

En este ejemplo, calculamos las proporciones de todas las combinaciones de los niveles de los factores, en lugar de cada nivel por separado.

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

wght_ppccombined<-function(sdata,cppc,wght,brr) {
nc<-1;
for (i in 1:length(cppc)) {
nc <- nc * length(levels(as.factor(sdata[,cppc[i]])));
}
mppc<-matrix(ncol=nc,nrow=2);
mppc[,]<-0;
cn<-c();
#####################################################################
# Para evitar calcular las combinaciones con un proceso recursivo, lo que
# requeriría otra función, usaremos el vector ccom para ir definiendo las
# combinaciones. Cada posición del vector se corresponderá con el factor
# correspondiente del vector cppc, y contendrá el ordinal de uno de los
# niveles del factor. En cada vuelta del bucle, el vector ccom definirá una
# combinación diferente de los niveles de todos los factores en cppc
######################################################################
ccom<-rep(1,length(cppc));
bw<-TRUE;
while(bw) {
cnc<-paste(names(sdata)[cppc[1]],
levels(as.factor(sdata[,cppc[1]]))[ccom[1]], sep="-");
for (i in 2:length(ccom)) {
cnc<-paste(cnc,names(sdata)[cppc[i]],
levels(as.factor(sdata[,cppc[i]]))[ccom[i]], sep="-");
}
cn<-c(cn, cnc);
for (i in length(cppc):1) {
if (ccom[i] < length(levels(as.factor(sdata[,cppc[i]])))) {
ccom[i] = ccom[i] + 1;
break;
}
else {
if (i == 1) {
bw<-FALSE;
break;
}
else {
ccom[i] = 1;
}
}
}
}
colnames(mppc)<-cn;
rownames(mppc)<-c("PPC","SE");
swght<-sum(sdata[,wght]);
ix<-1;
######################################################################
# volvemos a realizar el procedimiento anterior para procesar las
# diferentes combinaciones de niveles de los factores
######################################################################
bw<-TRUE;
ccom<-rep(1,length(cppc));
while(bw) {
rfact<-sdata[,cppc[1]]==levels(as.factor(sdata[,cppc[1]]))[ccom[1]];
for (i in 2:length(ccom)) {
rfact<-rfact & (sdata[,cppc[i]] ==
levels(as.factor(sdata[,cppc[i]]))[ccom[i]]);
}
mppc[1,ix]<-sum(sdata[rfact,wght]) / swght;
for (k in 1:length(brr)) {
sbrr<-sum(sdata[,brr[k]]);
ppcbrr<-sum(sdata[rfact,brr[k]]) / sbrr;
mppc[2,ix]<-mppc[2,ix] + (ppcbrr-mppc[1,ix])^2;
}
ix<-ix + 1;
for (i in length(cppc):1) {
if (ccom[i] < length(levels(as.factor(sdata[,cppc[i]])))) {
ccom[i] = ccom[i] + 1;
break;
}
else {
if (i == 1) {
bw<-FALSE;
break;
}
else {
ccom[i] = 1;
}
}
}
}
mppc[2,]<-sqrt((mppc[2,] * 4) / length(brr));
return(mppc);
}

El resultado es una matriz con dos filas, con una columna por cada una de las combinaciones de niveles de todos los factores y una fila para las proporciones y otra para los errores estándar.

Ejemplo 6, cálculo de las proporciones del cruce de dos grupos de factores

En este ejemplo, cruzaremos los valores de los niveles de un conjunto de factores con los de otro conjunto diferente y haremos el recuento de los casos ponderado por sus pesos, junto con el cálculo de los errores estandar.

El nombre de la función es wght_ppccrossed, y este es el código:

wght_ppccrossed<-function(sdata,rppc,cppc,wght,brr) {
nc<-0;
for (i in 1:length(cppc)) {
nc <- nc + length(levels(as.factor(sdata[,cppc[i]])));
}
nr<-0;
for (i in 1:length(rppc)) {
nr <- nr + length(levels(as.factor(sdata[,rppc[i]])));
}
mppc<-matrix(ncol=nc,nrow=nr);
mppc[,]<-0;
seppc<-matrix(ncol=nc,nrow=nr);
seppc[,]<-0;
result<-vector('list',2);
names(result)<-c("PPC","SE");
cn<-c();
for (i in 1:length(cppc)) {
for (j in 1:length(levels(as.factor(sdata[,cppc[i]])))) {
cn<-c(cn, paste(names(sdata)[cppc[i]],
levels(as.factor(sdata[,cppc[i]]))[j], sep="-"));
}
}
colnames(mppc)<-cn;
colnames(seppc)<-cn;
cr<-c();
for (i in 1:length(rppc)) {
for (j in 1:length(levels(as.factor(sdata[,rppc[i]])))) {
cr<-c(cr, paste(names(sdata)[rppc[i]],
levels(as.factor(sdata[,rppc[i]]))[j], sep="-"));
}
}
rownames(mppc)<-cr;
rownames(seppc)<-cr;
ir<-1;
swght<-sum(sdata[,wght]);
for (r in 1:length(rppc)) {
for (l in 1:length(levels(as.factor(sdata[,rppc[r]])))) {
rfact <- (sdata[,rppc[r]] ==
levels(as.factor(sdata[,rppc[r]]))[l]);
ic<-1;
for (i in 1:length(cppc)) {
for (j in 1:length(levels(as.factor(sdata[,cppc[i]])))) {
cfact<-(sdata[,cppc[i]] ==
levels(as.factor(sdata[,cppc[i]]))[j]) & rfact;
mppc[ir,ic]<-sum(sdata[cfact,wght]) / swght;
for (k in 1:length(brr)) {
sbrr<-sum(sdata[,brr[k]]);
ppcbrr<-sum(sdata[cfact,brr[k]]) / sbrr;
seppc[ir,ic]<-seppc[ir,ic] + (ppcbrr-mppc[ir,ic])^2;
}
ic<-ic + 1;
}
}
ir<-ir + 1;
}
}
seppc[,]<-sqrt((seppc[,] * 4) / length(brr));
result[[1]]<-mppc;
result[[2]]<-seppc;
return(result);
}

En el parámetro cppc pasaremos un vector con los índices o los nombres de las columnas de sdata que queremos cruzar con las columnas indicadas en el parámetro rppc, también con un vector de índices o nombres de columna.

Obtendremos una lista con dos matrices, la primera con las proporciones, entre 0 y 1, del cruce entre las filas (en rppc) y las columnas (en cppc) de los diferentes valores de los niveles de los factores, la segunda con los errores estándar correspondientes.

Ejemplo 7, cálculo simple de las proporciones de una serie de factores en una muestra con varios países

Este ejemplo es como el ejemplo 4, pero esta vez tenemos una muestra con más de un país.

El nombre de la función es wght_ppc_bycnt, y este es el código:

wght_ppc_bycnt<-function(sdata,cnt,cppc,wght,brr) {
nc<-0;
for (i in 1:length(cppc)) {
nc <- nc + length(levels(as.factor(sdata[,cppc[i]])));
}
nr<-length(levels(as.factor(sdata[,cnt])));
mppc<-matrix(ncol=nc,nrow=nr);
mppc[,]<-0;
seppc<-matrix(ncol=nc,nrow=nr);
seppc[,]<-0;
result<-vector('list',2);
names(result)<-c("PPC","SE");
cn<-c();
for (i in 1:length(cppc)) {
for (j in 1:length(levels(as.factor(sdata[,cppc[i]])))) {
cn<-c(cn, paste(names(sdata)[cppc[i]],
levels(as.factor(sdata[,cppc[i]]))[j], sep="-"));
}
}
colnames(mppc)<-cn;
colnames(seppc)<-cn;
cr<-c();
for (j in 1:length(levels(as.factor(sdata[,cnt])))) {
cr<-c(cr, levels(as.factor(sdata[,cnt]))[j]);
}
rownames(mppc)<-cr;
rownames(seppc)<-cr;
ir<-1;
for (l in 1:length(levels(as.factor(sdata[,cnt])))) {
rfact <- (sdata[,cnt]==levels(as.factor(sdata[,cnt]))[l]);
swght<-sum(sdata[rfact,wght]);
ic<-1;
for (i in 1:length(cppc)) {
for (j in 1:length(levels(as.factor(sdata[,cppc[i]])))) {
cfact<-(sdata[,cppc[i]] ==
levels(as.factor(sdata[,cppc[i]]))[j]) & rfact;
mppc[ir,ic]<-sum(sdata[cfact,wght]) / swght;
for (k in 1:length(brr)) {
sbrr<-sum(sdata[rfact,brr[k]]);
ppcbrr<-sum(sdata[cfact,brr[k]]) / sbrr;
seppc[ir,ic]<-seppc[ir,ic] + (ppcbrr-mppc[ir,ic])^2;
}
ic<-ic + 1;
}
}
ir<-ir + 1;
}
seppc[,]<-sqrt((seppc[,] * 4) / length(brr));
result[[1]]<-mppc;
result[[2]]<-seppc;
return(result);
}

El índice o nombre de la columna con el país se debe pasar en el parámetro cnt. El resultado es una lista con dos matrices. La primera tiene una fila por país y una columna con las proporciones de cada uno de los niveles de los factores, la segunda contiene los errores estándar correspondientes.

Ejemplos de cálculo de diferencias de medias y su error típico

En estos ejemplos podemos ver cómo realizar el cálculo de una diferencia entre dos medias y su error típico. En este enlace puedes descargar los ejemplos de código para calcular diferencias de medias.

En estos ejemplos el peso del estudiante y los pesos replicados se pasan en los parámetros wght y brr, como en los ejemplos anteriores.

Ejemplo 8, cálculo de la diferencia entre las medias de un conjunto de columnas agrupando por los niveles de un conjunto de factores

En este ejemplo, calcularemos la diferencia entre las medias para cada una de las combinaciones de dos de los niveles de un conjunto de factores. Un factor puede ser por ejemplo el sexo, y se calculará la diferencia entre las medias de un determinado índice entre los alumnos y las alumnas. Otros factores pueden tener más de dos niveles, y se calculará la diferencia de medias para cada combinación de dos de dichos niveles.

La función es wght_meandiff, y su código es:

wght_meandiff<-function(sdata,cmean,cfact,wght,brr) {
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;
}
}
}
mmeans<-matrix(ncol=nc,nrow=2 * length(cmean));
mmeans[,]<-0;
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="-"));
}
}
}
colnames(mmeans)<-cn;
rn<-c();
for (i in 1:length(cmean)) {
rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF",sep="-"));
rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF","SE",sep="-"));
}
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]])))) {
ir<-1;
rfact1<-sdata[,cfact[f]] ==
levels(as.factor(sdata[,cfact[f]]))[l];
rfact2<-sdata[,cfact[f]] ==
levels(as.factor(sdata[,cfact[f]]))[k];
swght1<-sum(sdata[rfact1,wght]);
swght2<-sum(sdata[rfact2,wght]);
for (i in 1:length(cmean)) {
mmeans[ir,ic]<-(sum(sdata[rfact1,wght] *
sdata[rfact1,cmean[i]])/swght1) -
(sum(sdata[rfact2,wght] *
sdata[rfact2,cmean[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,cmean[i]])/sbrr1) -
(sum(sdata[rfact2,brr[j]] *
sdata[rfact2,cmean[i]])/sbrr2);
mmeans[ir+1,ic]<-mmeans[ir+1,ic] + (mmbrj -
mmeans[ir,ic])^2;
}
mmeans[ir+1,ic]<-sqrt((mmeans[ir+1,ic] * 4) / length(brr));
ir<-ir+2;
}
ic<-ic + 1;
}
}
}
return(mmeans);
}

En cmean pasaremos un vector con los índices o los nombres de las columnas de las que queremos calcular la media, mientras que en cfact pasaremos un vector de índices o nombres de columnas con factores para agrupar los datos en función de sus niveles.

El resultado será una matriz con una columna para la diferencia de medias para cada combinación de dos niveles de cada uno de los factores, y dos filas para cada uno de los campos de los que queremos calcular la media, la primera para las diferencias y la segunda para los errores típicos de las mismas.

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

En este ejemplo veremos un ejemplo de cálculo de la diferencia de medias de una serie de columnas entre combinaciones de dos países de entre todos los que componen la muestra.

El nombre de la función es wght_meandiffbycnt1 y el código es el siguiente:

wght_meandiffbycnt1<-function(sdata,cmean,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 * length(cmean));
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(names(sdata)[cnt],
levels(as.factor(sdata[,cnt]))[j],
levels(as.factor(sdata[,cnt]))[k], sep="-"));
}
}
colnames(mmeans)<-cn;
rn<-c();
for (i in 1:length(cmean)) {
rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF",sep="-"));
rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF","SE",sep="-"));
}
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])))) {
ir<-1;
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]);
for (i in 1:length(cmean)) {
mcnt1<-sum(sdata[rcnt1,wght]*sdata[rcnt1,cmean[i]])/swght1;
mcnt2<-sum(sdata[rcnt2,wght]*sdata[rcnt2,cmean[i]])/swght2;
mmeans[ir,ic]<-mcnt1 - mcnt2;
secnt1<-0;
secnt2<-0;
for (j in 1:length(brr)) {
sbrr1<-sum(sdata[rcnt1,brr[j]]);
sbrr2<-sum(sdata[rcnt2,brr[j]]);
mbcnt1<-sum(sdata[rcnt1,brr[j]] *
sdata[rcnt1,cmean[i]])/sbrr1;
mbcnt2<-sum(sdata[rcnt2,brr[j]] *
sdata[rcnt2,cmean[i]])/sbrr2;
secnt1<-secnt1 + (mcnt1 - mbcnt1)^2;
secnt2<-secnt2 + (mcnt2 - mbcnt2)^2;
}
secnt1<-sqrt((secnt1 * 4) / length(brr));
secnt2<-sqrt((secnt2 * 4) / length(brr));
mmeans[ir+1,ic]<-sqrt(secnt1^2 + secnt2^2);
ir<-ir+2;
}
ic<-ic + 1;
}
}
return(mmeans);
}

En el parámetro cnt pasaremos el índice o el nombre de la columna con el país, y en cmean un vector con los índices o los nombres de las columnas para las que deseamos obtener las diferencias de medias. Lo que obtenemos es una matriz con una columna para la diferencia de medias entre cada combinación de dos países y dos filas para cada una de las columnas indicadas en cmean, la primera para la diferencia de medias, y la segunda para el error estándar de la diferencia.

Cuando se calcula el error estándar de una diferencia entre dos países, el cálculo es diferente al realizado con datos pertenecientes a un mismo país. En este caso, las submuestras de los dos países son independientes, y el error típico de la diferencia es igual a la raíz cuadrada de la suma de los errores típicos al cuadrado de cada una de las medias. Así se ha calculado el error típico en este ejemplo.

Ejemplo 10, cálculo de las diferencias entre medias de un conjunto de países, agrupadas por los niveles de un conjunto de factores

En este último ejemplo, también tenemos una muestra de varios países, pero lo que vamos a calcular son las diferencias de medias entre una serie de columnas según la combinación de niveles de una serie de factores tomados dos a dos, como en el ejemplo 4.

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

wght_meandiffbycnt2<-function(sdata,cmean,cnt,cfact,wght,brr) {
lcntrs<-vector('list',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];
}
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();
for (i in 1:length(cmean)) {
rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF",sep="-"));
rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF","SE",sep="-"));
}
for (p in 1:length(levels(as.factor(sdata[,cnt])))) {
mmeans<-matrix(ncol=nc,nrow=2 * length(cmean));
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]])))) {
ir<-1;
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]);
for (i in 1:length(cmean)) {
mmeans[ir,ic]<-(sum(sdata[rfact1,wght] *
sdata[rfact1,cmean[i]])/swght1) -
(sum(sdata[rfact2,wght] *
sdata[rfact2,cmean[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,cmean[i]])/sbrr1) -
(sum(sdata[rfact2,brr[j]] *
sdata[rfact2,cmean[i]])/sbrr2);
mmeans[ir+1,ic]<-mmeans[ir+1,ic] + (mmbrj -
mmeans[ir,ic])^2;
}
mmeans[ir+1,ic]<-sqrt((mmeans[ir+1,ic] * 4) /
length(brr));
ir<-ir+2;
}
ic<-ic + 1;
}
}
}
lcntrs[[p]]<-mmeans;
}
return(lcntrs);
}

Este ejemplo nos devuelve una lista con una matriz por país. La estructura de estas matrices es la misma que la del ejemplo 4.