Example R programs and commands 19. Chi-squared and log-likelihood goodness of fit tests # All lines preceded by the "#" character are my comments. # All other left-justified lines are my input. # All other indented lines are the R program output. # CHI-SQUARED test syntax # # chisq.test(x, y = NULL, correct = TRUE, # p = rep(1/length(x), length(x)), rescale.p = FALSE, ...) # # Arguments # x a vector or matrix. # y a vector; ignored if x is a matrix. # correct TRUE ==> apply continuity correction for 2x2 tables. # p vector of probabilities, same length as x. # rescale.p TRUE ==> p is rescaled (if necessary) to sum to 1. # FALSE ==> if p does not sum to 1, then complain. ####################################################################### # # Goodness-of-fit: # # Test if an observed frequency table is samples from a uniform density: # obs1 <- c(29, 23, 17, 25, 22, 15, 27) # observed frequencies chisq.test(obs1) # default p=c(1/7,...,1/7) is appropriate here Chi-squared test for given probabilities data: obs X-squared = 6.8987, df = 6, p-value = 0.3303 # ==> do not reject H0 at the .05 level. # Test if another observed frequency table is samples from a uniform density: obs2 <- c(29, 23, 10, 25, 22, 9, 27) # observed frequencies chisq.test(obs2) # default p=c(1/7,...,1/7) is appropriate here Chi-squared test for given probabilities data: obs X-squared = 18.6069, df = 6, p-value = 0.004882 # ==> reject H0 at the .05 level. ####################################################################### # # LOG-LIKELIHOOD computation # # Extract the observed and expected values. obs2 <- c(29, 23, 10, 25, 22, 9, 27) # observed frequencies fit<-chisq.test(obs2) # default p=c(1/7,...,1/7) is appropriate here summary(fit) Length Class Mode statistic 1 -none- numeric # Chi squared (X2) parameter 1 -none- numeric # Degrees of freedom (df) p.value 1 -none- numeric method 1 -none- character data.name 1 -none- character observed 7 -none- numeric expected 7 -none- numeric residuals 7 -none- numeric fit\$expected [1] 20.71429 20.71429 20.71429 20.71429 20.71429 20.71429 20.71429 fit\$observed [1] 29 23 10 25 22 9 27 fit\$observed/fit\$expected [1] 1.4000000 1.1103448 0.4827586 1.2068966 1.0620690 0.4344828 1.3034483 log(fit\$observed/fit\$expected) [1] 0.33647224 0.10467062 -0.72823850 0.18805223 0.06021886 -0.83359902 [7] 0.26501327 g<-2*sum(fit\$observed*log(fit\$observed/fit\$expected)) g [1] 21.12364 # Compute the p-value of g assuming it has the Chi-squared density: nu<-6; 1-pchisq(g,nu) [1] 0.001743084 ####################################################################### # # Heterogeneity Chi Squared (Method 1): # # Test if two frequency tables are samples from the uniform density # obs1 <- c(29, 23, 17, 25, 22, 15, 27) # observed frequencies, table 1 obs2 <- c(29, 23, 10, 25, 22, 9, 27) # observed frequencies, table 2 # # Individual chi-squareds sum to give the total chi-squared deviation # from the expected density p={1/K,...,1/K}: fit1<-chisq.test(obs1) fit2<-chisq.test(obs2) total.chisq <- fit1\$statistic + fit2\$statistic total.df <- fit1\$parameter + fit2\$parameter # Pooled chi-squared deviation from expected density p={1/K,...,1/K}: fitp<-chisq.test(obs1+obs2) pooled.chisq <- fitp\$statistic; pooled.chisq pooled.df <- fitp\$parameter; pooled.df # Compare deviations from expected: heterogeneity.chisq <- total.chisq - pooled.chisq; heterogeneity.chisq heterogeneity.df <- total.df - pooled.df; heterogeneity.df pchisq( heterogeneity.chisq, heterogeneity.df, lower.tail=FALSE) # NOTE: the uniform density "p=c(1/K,...,1/K)" with K=length(x) is R's # default, but any expected density may be used by specifying # "p=c(p1,p2,...pK)" in chisq.test(). ####################################################################### # # Heterogeneity Chi Squared (Method 2 -- Pearson's Chi Squared): # # Test if two frequency tables are samples from the pooled density # obs1 <- c(29, 23, 17, 25, 22, 15, 27) # observed frequencies, table 1 obs2 <- c(29, 23, 10, 25, 22, 9, 27) # observed frequencies, table 2 # # Test if the frequency tables are linearly dependent: fit <- chisq.test(cbind(obs1,obs2)); fit # NOTE: this is the same as the Method 1 test with the expected # density p=c(p1,...,pK) determined by the pooled frequencies. pooled.density <- (obs1+obs2)/sum(obs1+obs2) fit1 <- chisq.test(obs1,p=pooled.density); fit1 fit2 <- chisq.test(obs2,p=pooled.density); fit2 fitp <- chisq.test(obs1+obs2,p=pooled.density); fitp # NOTE: the pooled chi squared is exactly 0, so there is no need to # subtract it: heterogeneity.chisq <- fit1\$statistic + fit2\$statistic; heterogeneity.chisq # ....but the pooled DF is nonzero and must still be subtracted: heterogeneity.df <- fit1\$parameter + fit2\$parameter - fitp\$parameter; heterogeneity.df # Find the p-value of this test statistic manually: pchisq( heterogeneity.chisq, heterogeneity.df, lower.tail=FALSE ) ####################################################################### # # Kolmogorov-Smirnov Goodness-of-Fit: # # Test if two frequency tables comes from the same density # freq1 <- c(22, 25, 27, 15, 26, 24, 19) # observed frequencies, table 1 freq2 <- c(26, 23, 10, 25, 18, 29, 27) # observed frequencies, table 2 # For R's K-S test, the inputs are observations not frequencies: x1<-rep(1:length(freq1),freq1); x1 # generate observations from freq1 x2<-rep(1:length(freq2),freq2); x2 # generate observations from freq2 ks.test(x1,x2) # test whether x1,x2 come from the same density # # Test if a sequence of observations comes from a specified density # # Specify the uniform density on the discrete set {1,2,...,7}: ks.test(x1,"punif",min=0.5,max=7.5) # {1,...,7} <===> [0.5,7.5] ks.test(x2,"punif",min=0.5,max=7.5) # {1,...,7} <===> [0.5,7.5] # # REMARK: if p is large so H_0 is not rejected, then the frequency # tables may be added (pooled.freq <- freq1+freq2) and the observations # may be concatenated (pooled.x <- c(x1,x2) ).