0
$\begingroup$

I have data from an experiment testing the number of 'cases' at each of three measurement points (0, 12, and 24 weeks). I am interested in whether the proportion of cases across these measurement periods is the same in two experimental groups (placebo vs treatment). Following @gung's excellent response here ($\chi^2 $ of multidimensional data) I created a 3d array.

counts <- c(37, 14, 10, 30, 26, 20, 31,  9,  8, 30, 28, 18)
threeDArray <- array(counts, dim = c(3,2,2))
tab <- as.table(threeDArray)
names(dimnames(tab)) <- c("week", "caseness", "treatment")
dimnames(tab)[[1]] = c("0", "12", "24")
dimnames(tab)[[2]] = c("case", "nonCase")
dimnames(tab)[[3]] = c("placebo", "treatment") 
tab

# , , treatment = placebo
# 
# caseness
# week case nonCase
# 0    37      30
# 12   14      26
# 24   10      20
# 
# , , treatment = treatment
# 
# caseness
# week case nonCase
# 0    31      30
# 12    9      28
# 24    8      18

In order to test whether the distribution of cases across weeks is the same in both experimental groups I created this model, where there is no treatment factor...

library(MASS)

m1 <- loglm(~week*caseness, tab)

# Statistics:
#                       X^2 df  P(> X^2)
# Likelihood Ratio 2.027869  6 0.9171177
# Pearson          2.017928  6 0.9180424

...and then this model where the treatment effect is added

m2 <- loglm(~treatment + week*caseness, tab)

# Statistics:
#                       X^2 df  P(> X^2)
# Likelihood Ratio 1.380092  5 0.9264786
# Pearson          1.373826  5 0.9271548

Apparently this..

summary(tab)
# Number of cases in table: 261 
# Number of factors: 3 
# Test for independence of all factors:
#   Chisq = 14.658, df = 7, p-value = 0.04064

...is the null model, which is that all cells have the same expected value, and this...

msat <- loglm(~week*caseness*treatment, tab)

# Statistics:
#                  X^2 df P(> X^2)
# Likelihood Ratio   0  0        1
# Pearson            0  0        1

Now back to models 1 and 2. When I run the likelihood ratio test...

anova(m1, m2)

#           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
# Model 1   2.027869  6                                    
# Model 2   1.380092  5  0.6477776         1        0.42091
# Saturated 0.000000  0  1.3800919         5        0.92648

...it compares m1 to m2 and then m2 to the saturated model.

From this output it looks like there no significant difference between m1 and m2 and no significant difference between either model and the saturated model

pchisq(2.027869, 6, lower.tail = F) # test difference in deviance between m1 and saturated

# [1] 0.9171177

Question 1: Why does m1 have 6 df and m2 only 5, when m2 has an extra variable?

Question 2: What is the difference between the null model and the saturated model? I assumed a model where the expected cell counts are all the same, and a saturated model, where all three factors are crossed with one another, were the same thing, but one (the null model) has 7 df and one has 0.

Question 3: How do I interpret these results? It looks to me, based on the likelihood ratio test, as if the inclusion of the treatment factor gives no real explanatory advantage and could be dropped from the model, but if neither m1 or m2 offer any significant improvement over the saturated model, why would one choose any model over any other?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.