name - title image matlab

Prueba de independencia de Matlab

Question

Para 1,000,000 de observaciones, observé un evento discreto, X, 3 veces para el grupo de control y 10 veces para el grupo de prueba.

Necesito realizar una prueba de independencia de Chi cuadrado en Matlab. Así es como lo harías en r:

m <- rbind(c(3, 1000000-3), c(10, 1000000-10))
#      [,1]   [,2] 
# [1,]    3 999997
# [2,]   10 999990
chisq.test(m)

La función r devuelve chi-squared = 2.7692, df = 1, p-value = 0.0961.

¿Qué función de Matlab debo usar o crear para hacer esto?

Answer 1

Aquí está mi propia implementación que uso:

function [hNull pValue X2] = ChiSquareTest(o, alpha)
    %#  CHISQUARETEST  Pearson''s Chi-Square test of independence
    %#
    %#    @param o          The Contignecy Table of the joint frequencies
    %#                      of the two events (attributes)
    %#    @param alpha      Significance level for the test
    %#    @return hNull     hNull = 1: null hypothesis accepted (independent)
    %#                      hNull = 0: null hypothesis rejected (dependent)
    %#    @return pValue    The p-value of the test (the prob of obtaining
    %#                      the observed frequencies under hNull)
    %#    @return X2        The value for the chi square statistic
    %#

    %# o:     observed frequency
    %# e:     expected frequency
    %# dof:   degree of freedom

    [r c] = size(o);
    dof = (r-1)*(c-1);

    %# e = (count(A=ai)*count(B=bi)) / N
    e = sum(o,2)*sum(o,1) / sum(o(:));

    %# [ sum_r [ sum_c ((o_ij-e_ij)^2/e_ij) ] ]
    X2 = sum(sum( (o-e).^2 ./ e ));

    %# p-value needed to reject hNull at the significance level with dof
    pValue = 1 - chi2cdf(X2, dof);
    hNull = (pValue > alpha);

    %# X2 value needed to reject hNull at the significance level with dof
    %#X2table = chi2inv(1-alpha, dof);
    %#hNull = (X2table > X2);

end

Y un ejemplo para ilustrar:

t = [3 999997 ; 10 999990]
[hNull pVal X2] = ChiSquareTest(t, 0.05)

hNull =
     1
pVal =
     0.052203
X2 =
       3.7693

Tenga en cuenta que los resultados son diferentes a los suyos porque chisq.test realiza una corrección por defecto, de acuerdo con ?chisq.test

correcto: un indicador lógico que indica si aplicar la corrección de continuidad al calcular el estadístico de prueba para tablas 2x2: una mitad se resta de todos | O - E | diferencias

Alternativamente, si tiene las observaciones reales de los dos eventos en cuestión, puede usar la función CROSSTAB que calcula la tabla de contingencia y devuelve las medidas de Chi2 y p-value:

X = randi([1 2],[1000 1]);
Y = randi([1 2],[1000 1]);
[t X2 pVal] = crosstab(X,Y)

t =
   229   247
   257   267
X2 =
     0.087581
pVal =
      0.76728

el equivalente en R sería:

chisq.test(X, Y, correct = FALSE)

Nota: Ambos enfoques (MATLAB) anteriores requieren la Caja de herramientas de estadísticas

Answer 2

Esta función probará la independencia utilizando la estadística de chi-cuadrado de Pearson y la estadística de razón de verosimilitud, junto con el cálculo de los residuos. Sé que esto se puede vectorizar más, pero estoy tratando de mostrar las matemáticas para cada paso.

function independenceTest(data)
df = (size(data,1)-1)*(size(data,2)-1); % Mean Degrees of Freedom
sd = sqrt(2*df);                        % Standard Deviation

u         = nan(size(data)); % Estimated expected frequencies 
p         = nan(size(data)); % Values used to calculate chi-square
lr        = nan(size(data)); % Values used to calculate likelihood-ratio
residuals = nan(size(data)); % Residuals

rowTotals    = sum(data,1);
colTotals    = sum(data,2);
overallTotal = sum(rowTotals);

%% Calculate estimated expected frequencies
for r=1:1:size(data,1)
    for c=1:1:size(data,2)
        u(r,c) = (rowTotals(c) * colTotals(r)) / overallTotal;
    end
end

%% Calculate chi-squared statistic
for r=1:1:size(data,1)
    for c=1:1:size(data,2)
        p(r,c) = (data(r,c) - u(r,c))^2 / u(r,c);
    end
end
chi = sum(sum(p)); % Chi-square statistic

%% Calculate likelihood-ratio statistic
for r=1:1:size(data,1)
    for c=1:1:size(data,2)
        lr(r,c) = data(r,c) * log(data(r,c) / u(r,c));
    end
end
G = 2 * sum(sum(lr)); % Likelihood-Ratio statisitc

%% Calculate residuals
for r=1:1:size(data,1)
    for c=1:1:size(data,2)
        numerator   = data(r,c) - u(r,c);
        denominator = sqrt(u(r,c) * (1 - colTotals(r)/overallTotal) * (1 - rowTotals(c)/overallTotal));
        residuals(r,c) = numerator / denominator;
    end
end