In Tuesday’s class we started talking about Likelihood Ratio Tests. In this test, we compare two nested models in terms of model fit. A nested model means that one model is a special case of another; for example, \({y} = {\beta_0}\) is a nested model of \({y} = {\beta_0} + {\beta_1}x\).
To get our test statistic for thsi test, we take the likelihood of one model and divide it by the likelihood of the other, and then take the log of it. In this test statistic, the likelihood of the complicated model should be larger than the likelihood of the simpler model.
The distribution of this test stat is Chi Squared. The degrees of freedom we use in this test is the difference in the number of parameters of the two models.
To obtain our p-value, we take area to the right of our test statistic because the p-value in this distribution is the probability of getting test statistic as extreme or more extreme than ours.
The steps for this test are:
If we get a test stat close to 0, this means that our complicated and simpler model are pretty similar, so we should stick with the simpler model.
As an example, we will use the iris dataset. We create two models with sepal length as the response, but one model will just have petal length as the predictor and the other will have petal length and species.
data(iris)
attach(iris)
mymod <- lm(Sepal.Length ~ Petal.Length)
mymod2 <- lm(Sepal.Length ~ Petal.Length + Species)
summary(mymod2)
##
## Call:
## lm(formula = Sepal.Length ~ Petal.Length + Species)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75310 -0.23142 -0.00081 0.23085 1.03100
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.68353 0.10610 34.719 < 2e-16 ***
## Petal.Length 0.90456 0.06479 13.962 < 2e-16 ***
## Speciesversicolor -1.60097 0.19347 -8.275 7.37e-14 ***
## Speciesvirginica -2.11767 0.27346 -7.744 1.48e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.338 on 146 degrees of freedom
## Multiple R-squared: 0.8367, Adjusted R-squared: 0.8334
## F-statistic: 249.4 on 3 and 146 DF, p-value: < 2.2e-16
anova(mymod, mymod2)
## Analysis of Variance Table
##
## Model 1: Sepal.Length ~ Petal.Length
## Model 2: Sepal.Length ~ Petal.Length + Species
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 148 24.525
## 2 146 16.682 2 7.8434 34.323 6.053e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Now that we’ve looked a bit at the comparison of the two models with the anova function, we obtain the log likelihood of each model using the logLik() function:
logLik(mymod)
## 'log Lik.' -77.02021 (df=3)
logLik(mymod2)
## 'log Lik.' -48.11637 (df=5)
We can see that the log likelihood for our simpler model is about -77.02 and the log likelihood of the more complex model is about -48.17.
We can now obtain our test stat to run our test. For this test, our null hypothesis is that we don’t need species in our model. Our alternative hypothesis is that we do need species in our model.
In this case, there are 2 parameters in our simpler model and 4 in our complex, so we have 2 degrees of freedom.
We obtain our test stat with the following equation:
teststat <- -2*(as.numeric(logLik(mymod))-as.numeric(logLik(mymod2)))
teststat
## [1] 57.80768
We now obtain our p-value:
pchisq(teststat, df = 2, lower.tail = FALSE)
## [1] 2.800409e-13
Because we get such a small p-value, we reject the null in favor of the alternative. We have really strong evidence that even after having petal length in the model, we should still have species in the model.
We will now run the same type of test on the time series models we created in Learning Log 23. We follow all the same steps as above but with the time series models:
library(forecast)
## Warning: package 'forecast' was built under R version 3.4.4
data(wineind)
mod1 <- arima(wineind, order = c(0,0,0), season = list( order = c( 1,1,1), period=12))
mod2 <- arima(wineind, order = c(0,0,0), season = list( order = c( 0,1,1), period=12))
logLik(mod1)
## 'log Lik.' -1520.259 (df=3)
logLik(mod2)
## 'log Lik.' -1520.341 (df=2)
teststat <- 2*(as.numeric(logLik(mod1))-as.numeric(logLik(mod2)))
teststat
## [1] 0.1644043
pchisq(teststat, df = 1, lower.tail = FALSE)
## [1] 0.6851336
Because we have a very high p-value, our test is telling us to stick with a simpler model. That is, we should stick with our model where P = 0, D = 1, Q = 1, versus the model where P = 1, D = 1, Q =1.
If I were to interpret these results to the wine salesman, I would tell them that this model definitely has a seasonal component. We can use our model to predict into the future:
par(mfrow=c(2,1))
forecast(mod2, h = 12)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 1994 23614.94 20335.48 26894.39 18599.448 28630.42
## Oct 1994 27679.00 24399.54 30958.45 22663.509 32694.48
## Nov 1994 32225.64 28946.19 35505.09 27210.152 37241.13
## Dec 1994 37271.66 33992.21 40551.11 32256.171 42287.15
## Jan 1995 14803.95 11524.49 18083.40 9788.458 19819.43
## Feb 1995 21919.62 18640.17 25199.07 16904.134 26935.11
## Mar 1995 23741.35 20461.90 27020.81 18725.868 28756.84
## Apr 1995 26136.40 22856.95 29415.85 21120.910 31151.88
## May 1995 24132.74 20853.28 27412.19 19117.248 29148.22
## Jun 1995 26534.20 23254.75 29813.65 21518.712 31549.69
## Jul 1995 29636.48 26357.02 32915.93 24620.989 34651.96
## Aug 1995 25267.08 21987.63 28546.53 20251.594 30282.57
Looking at this prediction table, we could give a wine salesman some advice about the future. For example, for August of 1995 (a year beyond this dataset), our suggestion would be to stock enough wine for about 25267 sales. To be safe, we would say that we are 80% confident that the total wine sales would be between 21987.63 and 28546.53, so the winesalesman would potentially need to stock even more than our point prediction.