First a quick review of an earlier topic, testing the hypothesis that the proportion of items in a population is equal to some value against an alternative hypothesis. (This was done in Topic 17.)

Let us say that we have a die which we suspect has been weighted so that the "4" comes up more often than it should. We know that in a fair die the probability of getting a "4" should be ⅙ , that is, in a huge population of rolls of that die the proportion of times that a "4" shows up should be really close to ⅙ . If we want to run a test of the null hypothesis

Let us do that. We will roll the die 96 times.

Now we need to find the number of 4's in the table. We can generate and verify the table in R and then count the number of 4's in R using

source("../gnrnd4.R") gnrnd4( key1=347129507, key2=6585376 ) L1 # now find the number of "4"s hold_table <- table( L1 ) hold_table # we are just interested in the number of 4's num_4s <- as.integer( hold_table[4] ) num_4sThis will produce the console output shown in Figure 1.

Knowing that "4" appears 23 times in our sample of 96 tosses means that we can test, at the 0.025 level of significance, the null hypothesis that the proportion of 4's is ⅙,

`hypoth_test_prop`

function,
source("../hypo_prop.R") hypoth_test_prop( 1/6, num_4s, 96, 7, 0.025)to produce

The result is that the attained significance is 0.0276 which is not less than our level of significance, 0.025. Therefore, we do not have sufficient evidence to reject

But, what about the other possible values, the "1", "2", "3", "5", and "6"? Should we run the same test, possibly changing "weighted in favor of" to "weighted against", for each of these values? We certainly have the ability to do this. It would just be running the

`hypoth_test_prop()`

function with the appropriate values for each alternative.
This leaves us with the following question:

First calculate all of this in R.

# get all of the observed values obs <- as.integer( hold_table ) obs sum(obs) # set the proportions to 1/6 props <- rep( 1/6, 6) props sum(props) # find the expected value for 96 tosses expected <- 96*props expected sum(expected)The console output for these is shown in Figure 3.

Table 1 | |||||||

item | 1 | 2 | 3 | 4 | 5 | 6 | total |

observed count | 25 | 7 | 10 | 23 | 15 | 16 | 96 |

H_{0} proportion | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1 |

expected | 16 | 16 | 16 | 16 | 16 | 16 | 96 |

# find the difference, i.e., the observed - expected values diff <- obs - expected diff sum( diff )

Table 2 | |||||||

item | 1 | 2 | 3 | 4 | 5 | 6 | total |

observed count | 25 | 7 | 10 | 23 | 15 | 16 | 96 |

H_{0} proportion | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1 |

expected | 16 | 16 | 16 | 16 | 16 | 16 | 96 |

observed – expected | 9 | -9 | -6 | 7 | -1 | 0 | 0 |

`observed - expected`

is interesting, it does not
tell us much, or at least, the sum of those values does not tell us much since it is always
going to be 0. Way back when we were learning about #find the square of the differences and the sum of those squares diff_sqr <- diff^2 diff_sqr sum( diff_sqr )

Table 3 | |||||||

item | 1 | 2 | 3 | 4 | 5 | 6 | total |

observed count | 25 | 7 | 10 | 23 | 15 | 16 | 96 |

H_{0} proportion | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1 |

expected | 16 | 16 | 16 | 16 | 16 | 16 | 96 |

observed – expected | 9 | -9 | -6 | 7 | -1 | 0 | 0 |

(observed – expected)² | 81 | 81 | 36 | 49 | 1 | 0 | 248 |

# now divide those values by the expected value component <- diff_sqr / expected component sum( component )

Table 4 | |||||||

item | 1 | 2 | 3 | 4 | 5 | 6 | total |

observed count | 25 | 7 | 10 | 23 | 15 | 16 | 96 |

H_{0} proportion | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1 |

expected | 16 | 16 | 16 | 16 | 16 | 16 | 96 |

observed – expected | 9 | -9 | -6 | 7 | -1 | 0 | 0 |

(observed – expected)² | 81 | 81 | 36 | 49 | 1 | 0 | 248 |

(observed – expected)²expected | 5.0625 | 5.0625 | 2.25 | 3.0625 | 0.0625 | 0 | 15.5 |

# find the critical high value qchisq( 0.025, 5, lower.tail = FALSE)

Using the critical value approach, 15.5 is larger than 12.8325 so we would reject

We have also used the attained significance approach to hypothesis testing problems. Here we got our sum of the values

# use the attained significance approach pchisq( 15.5, 5, lower.tail=FALSE)

As always, if our attained significance level is smaller than our stated level of significance then we reject the null hypothesis. Indeed, 0.0084 is less than 0.025 so we would reject

Before we move on, we need to re-examine Table 4. That table is organized as rows where we have columns for each of the observed items, plus a column for the totals. We will transpose that table, giving exactly the same values, but organizing it to have rows for each of the observed items.

Table 4 (transposed) | ||||||

item | obs count |
H_{0} proportion |
expected | obs-expected | (obs-expected)² | (obs-expected)² expected |

1 | 25 | 1/6 | 16 | 9 | 81 | 5.0625 |

2 | 7 | 1/6 | 16 | -9 | 81 | 5.0625 |

3 | 10 | 1/6 | 16 | -6 | 36 | 2.2500 |

4 | 23 | 1/6 | 16 | 7 | 49 | 3.0625 |

5 | 15 | 1/6 | 16 | -1 | 1 | 0.0625 |

6 | 16 | 1/6 | 16 | 0 | 0 | 0.0000 |

TOTAL | 96 | 1 | 96 | 0 | 248 | 15.5000 |

# make a data frame to hold most of our table our_table <- data.frame( 1:6, obs, rep(1/6,6), expected, diff, diff_sqr, component) colnames( our_table ) <- c("Items","obs","prop","expected", "diff", "diff_sqr","component") our_table

Before we talk about restrictions to this methodology and before we try to automate this process, let us look at a new problem.

Here is a new problem: We have a funny shaped solid that has 9 flat sides. The sides are numbered from 1 to 9. We have a hypothesis that probability that the solid, when "tossed", shows a particular number is given by P(X=1)=5/53, P(X=2)=7/53, P(X=3)=3/53, P(X=4)=8/53, P(X=5)=4/53, P(X=6)=6/53, P(X=7)=9/53, P(X=8)=3/53, and P(X=9)=8/53. We would run this test, say at the 0.05 level, the same way that we ran the test above. First, set up the probabilities.

props <- c(5,7,3,8,4,6,9,3,8)/53Then get our sample.

source( "../gnrnd5.R") gnrnd5( 77362025307, 4386483869) head(L1,10) tail(L1,10)

Then find out how many times we observe each of the nine faces.

hold_table <- table( L1 ) hold_table # get all of the observed values obs <- as.integer( hold_table ) obs sum(obs)

Then find the expected values along with the differences, observed - expected.

# find the expected value for 254 tosses expected <- 254*props expected sum(expected) # find the difference, i.e., the observed - expected values diff <- obs - expected diff sum(diff )

Again, when we add up the differences, we get 0, or really close to it given the rounding errors in the computer.

Moving forward we find the squares of the differences and then we find the quotient of those squared values and the expected values, along with finding the sum of those quotients.

#find the square of the differences and the sum of those squares diff_sqr <- diff^2 diff_sqr sum( diff_sqr ) # now divide those values by the expected value component <- diff_sqr / expected component sum( component )

Noting that we have 9 different categories, we know that we will have 8 degrees of freedom. We can find our critical value for our 0.05 level of significance test via

# find the critical high value qchisq( 0.05, 8, lower.tail = FALSE)

Thus, our test statistic, 15.81586 is greater than our critical high value 15.50731 so we reject the null hypothesis at the 0.05 level of significance. Alternatively, we could have used the attained significance approach and found out how strange it would be, if the null hypothesis were true, to get a sample with results showing that sum of 15.81586?

# use the attained significance approach pchisq( 15.81586, 8, lower.tail=FALSE)

The attained significance is 0.045 which is less than our stated level of 0.05 so we would reject the null hypothesis and say that we have statistical evidence that the funny solid shape does not have the proportion of faces on tosses as that was given in the problem statement.

`goodfit()`

does this.
# Roger Palay copyright 2021-04-12 # Saline, MI 48176 # # in the goodfit function: # items is a list of the item values (i.e., labels) # H_null is a list of the hypothesized proportions # observed is a list of the observed counts of items # sig_level is the level at which the test is to be run # show_table if TRUE will cause the function to use # View to display a nice table of values goodfit <- function( items, H_null, observed, sig_level,show_table=FALSE ) { # there are a few checks that we should make # first, is the sum of the H_null proportions == 1 sum_H0 <- sum( H_null ) if( abs( sum_H0 - 1) > 0.00000001) { return("H_null proportions must add to 1")} # second, all of the H_null values must be non-negative find_negs <- which( H_null < 0 ) if( length( find_negs ) > 0 ) {return( "cannot have negative proportions")} # Third, we need the same number of items, # H_null proportions, and observed frequencies items_length <- length( items ) H0_length <- length( H_null ) obs_length <- length( observed ) if( items_length != H0_length | items_length != obs_length) { return("Need same number of items, proportions, and observed frequencies")} obs_total <- sum(observed) expected <- H_null*obs_total # Fourth, we need enough expected values num_too_small <- length( which( expected < 5)) if( num_too_small > 0 ) { cat( "expected values, some too small","\n") cat(expected,"\n") return("use larger sample")} diff <- observed - expected diff2 <- diff^2 chi_component <- diff2/expected total = sum( chi_component ) goff_table <<- data.frame( items, H_null, observed, expected, diff, diff2, chi_component) colnames( goff_table ) <- c( "Items","H_null", "observed", "expected", "diff", "diff^2", "chi component") print( goff_table) if( show_table ) {View( goff_table )} cat( "-----------------------------------\n") cat("total observations = ",obs_total,"\n") cat("Number of categoies = ", obs_length,"\n") cat("Number of degrees of freedom = ", (obs_length-1),"\n") cat("Total of chi component =", total,"\n") cat("P(total >= ",total,") = ", pchisq(total, obs_length-1,lower.tail=FALSE), "\n") crit_value <- qchisq( sig_level, obs_length-1, lower.tail=FALSE) cat("For ",(sig_level*100),"% to the right, the ","\n") cat(" critical value is ", crit_value, "\n") if( total > crit_value ) { "reject H0"} else {" not enough evidence to reject H0"} }You might notice that the first parameter that you need to give the function is a list of the names of the categories. For the numbers 1 to 9 this is as easy as using

hold_colors <- c("blue","red","green","orange","purple")and then use that variable as the first parameter of the function call.

The second parameter gives the proportions (probabilities) for the null hypothesis. In the previous problem we put those into a variable called

`props`

via
props <- c(5,7,3,8,4,6,9,3,8)/53Therefore, we could use

`props`

as the second parameter.
The third parameter is our list of the number of times each of the categories is observed in the sample. Above we had taken a few steps to make list:

hold_table <- table( L1 ) hold_table # get all of the observed values obs <- as.integer( hold_table ) obsIf we have done that then we could use

`obs`

as the third parameter.
However, if you were looking for a shorcut, you could just use `as.integer(table(L1))`

in that place.
The fourth parameter is just the level of significance for our test. And the fifth parameter, which is optional, specifies if you want the function to automatically force a

`goff_table`

, in the environment.
Even if you do not ask the function to show the data frame it is still there for you.
So, assuming that you have generated the sample in L1, then you could perform our last test via

source("../goodfit.r") props <- c(5,7,3,8,4,6,9,3,8)/53 goodfit(1:9, props, as.integer(table(L1)), 0.05,FALSE)And this produces the output shown in Figure 16.

Below are the R lines of script that were used to generate this page:

# statements used in the production of # the goodnessoffit.htm web page # # generate the data for Table 1 source("../gnrnd4.R") gnrnd4( key1=347129507, key2=6585376 ) L1 # now find the number of "4"s hold_table <- table( L1 ) hold_table # we are just interested in the number of 4's num_4s <- as.integer( hold_table[4] ) num_4s # we are ready to run our test source("../hypo_prop.R") hypoth_test_prop( 1/6, num_4s, 96, 7, 0.025) # get all of the observed values obs <- as.integer( hold_table ) obs sum(obs) # set the proportions to 1/6 props <- rep( 1/6, 6) props sum(props) # find the expected value for 96 tosses expected <- 96*props expected sum(expected) # find the difference, i.e., the observed - expected values diff <- obs - expected diff sum(diff ) #find the square of the differences and the sum of those squares diff_sqr <- diff^2 diff_sqr sum( diff_sqr ) # now divide those values by the expected value component <- diff_sqr / expected component sum( component ) # find the critical high value qchisq( 0.025, 5, lower.tail = FALSE) # use the attained significance approach pchisq( 15.5, 5, lower.tail=FALSE) # make a data frame to hold most of our table our_table <- data.frame( 1:6, obs, rep(1/6,6), expected, diff, diff_sqr, component) colnames( our_table ) <- c("Items","obs","prop","expected", "diff", "diff_sqr","component") our_table ####################################### # here is a new problem: We have a funny shaped # solid that has 9 flat sides. The sides are numbered # from 1 to 9. We have a hypothesis that probability # that the solid, when "tossed", shows a particular number is # given by P(X=1)=5/53, P(X=2)=7/53, P(X=3)=3/53, P(X=4)=8/53, # P(X=5)=4/53, P(X=6)=6/53, P(X=7)=9/53, P(X=8)=3/53, and # P(X=9)=8/53. props <- c(5,7,3,8,4,6,9,3,8)/53 # Let us get a sample of such tosses and test the hypothesis # that those proportions are correct against the alternative # that they are not. source( "../gnrnd5.R") gnrnd5( 77362025307, 4386483869) head(L1,10) tail(L1,10) hold_table <- table( L1 ) hold_table # get all of the observed values obs <- as.integer( hold_table ) obs sum(obs) # find the expected value for 254 tosses expected <- 254*props expected sum(expected) # find the difference, i.e., the observed - expected values diff <- obs - expected diff sum(diff ) #find the square of the differences and the sum of those squares diff_sqr <- diff^2 diff_sqr sum( diff_sqr ) # now divide those values by the expected value component <- diff_sqr / expected component sum( component ) # find the critical high value qchisq( 0.05, 8, lower.tail = FALSE) # use the attained significance approach pchisq( 15.81586, 8, lower.tail=FALSE) # now show off the function to do this source("../goodfit.r") props <- c(5,7,3,8,4,6,9,3,8)/53 goodfit(1:9, props, as.integer(table(L1)), 0.05,FALSE)

©Roger M. Palay Saline, MI 48176 April, 2020