Source: vignettes/inference.Rmd

Foehnix Demo

foehnix objects also provide asymptotic inference of the estimated coefficients of the mixture model (see statistical model and logistic regression with IWLS).

Let’s start with estimating our demo flexmix foehn diagnosis model. We have already prepared the data sets (object data). Details how to generate the data object and more information about the foehnix model specification can be found on the getting started manual page.

# Estimate the model
mod <- foehnix(diff_t ~ ff + rh, data = data)

The summary method prints the test statistics for both parts of the mixture model, the component model and the concomitant model (if specified).

summary(mod, detailed = TRUE)
## 
## Call: foehnix(formula = diff_t ~ ff + rh, data = data)
## 
## Number of observations (total)   113952 (5527 due to inflation)
## Removed due to missing values     20014 (17.6 percent)
## Outside defined wind sector           0 (0.0 percent)
## Used for classification           93938 (82.4 percent)
## 
## Climatological foehn occurance 75.87 percent (on n = 93938)
## Mean foehn probability 75.61 percent (on n = 93938)
## 
## Log-likelihood: -221461.2, 7 effective degrees of freedom
## Corresponding AIC = 442936.3, BIC = 443002.5
## Time required for model estimation: 19.8 seconds
## 
## Cluster separation (ratios close to 1 indicate
## well separated clusters):
##                        prior  size post>0 ratio
## Component 1 (foehn)     0.76 71268  93682  0.76
## Component 2 (no foehn)  0.24 22670  62085  0.37
## 
## ---------------------------------
## 
## Concomitant model: z test of coefficients
##                   Estimate  Std. error  z value  Pr(>|z|)    
## cc.(Intercept)  2.34150441  0.04768558   49.103 < 2.2e-16 ***
## cc.ff          -1.88353356  0.00519504 -362.564 < 2.2e-16 ***
## cc.rh           1.76129854  0.00077906 2260.790 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Number of IWLS iterations: 1 (algorithm converged)
## Dispersion parameter for binomial family taken to be 1.

Component Model Inference

The inference for the two location parameters μ1\mu_1 and μ2\mu_2 of the two components are based on the asymptotic theory. We expect that the estimates of our coefficients are unbiased. Thus, the expectation of our estimated coefficients is the estimated coefficient itself (E(μ̂1)=μ1\text{E}(\hat{\mu}_1) = \mu_1, E(μ̂2)=μ2\text{E}(\hat{\mu}_2) = \mu_2).

In a general form the covariance matrix of a liner model (one component, unweighted) for a set of i=1,…,Ni = 1, \dots, N observations can be expressed as follows:

  • Cov(β)=1N−P𝜖⊤𝜖(𝐗⊤𝐗)−1\text{Cov}(\beta) = \frac{1}{N - P} \mathit{\epsilon}^\top \mathit{\epsilon} (\mathbf{X}^\top \mathbf{X})^{-1}

… where NN is the sample size, PP the number of parameters or covariates, 𝜖=𝑦−μ•\mathit{\epsilon} = \mathit{y} - \mu_\bullet the model residuals, and 𝐗\mathbf{X} the model matrix of the linear model. In case of a foehnix model we have two components where each component consists of an intercept only model (μ1\mu_1 and μ2\mu_2 do not depend on additional covariates). Thus, P=1P = 1 and 𝐗\mathbf{X} is an N×1N \times 1 matrix with 1s. The estimates are based on a set of weighted observation yy, in this example diff_t (y=y =diff_t). The weights are the a-posteriori probabilities 𝑝̂\hat{\mathit{p}} (the foehn probabilities) of the foehnix model and have to be taken into account when calculating the standard errors. With these weights the standard error for component 2 can be written as:

  • sd̂(μ2)2=1∑i=1Np̂i−1(𝜖⋅𝑝̂)⊤(𝜖⋅𝑝̂)((𝐗⋅𝑝̂)⊤(𝐗⋅𝑝̂))−1\hat{\text{sd}}(\mu_2)^2 = \frac{1}{\sum_{i=1}^N \hat{p}_i - 1} (\mathit{\epsilon} \cdot \hat{\mathit{p}})^\top (\mathit{\epsilon} \cdot \hat{\mathit{p}}) \big((\mathbf{X} \cdot \hat{\mathit{p}})^\top (\mathbf{X} \cdot \hat{\mathit{p}})\big)^{-1}

… or much simpler:

  • sd̂(μ2)2=1∑i=1Np̂i2⋅(∑i=1Np̂i−1)∑i=1N((yi−μ2)⋅p̂i)2\hat{\text{sd}}(\mu_2)^2 = \frac{1}{\sum_{i=1}^N \hat{p}_i^2 \cdot (\sum_{i=1}^N \hat{p}_i - 1)} \sum_{i=1}^N \big((y_i - \mu_2) \cdot \hat{p}_i\big)^2

The same holds for component 1 except that our weights are 1−𝑝̂1 - \hat{\mathit{p}}:

  • sd̂(μ1)2=1∑i=1N(1−p̂i)2⋅(∑i=1N(1−p̂i)−1)∑i=1N((yi−μ2)⋅(1−p̂i))2\hat{\text{sd}}(\mu_1)^2 = \frac{1}{\sum_{i=1}^N (1 - \hat{p}_i)^2 \cdot (\sum_{i=1}^N (1 - \hat{p}_i) - 1)} \sum_{i=1}^N \big((y_i - \mu_2) \cdot (1 - \hat{p}_i)\big)^2

Concomitant Model Inference

If a concomitant model has been specified summary will also return the corresponding z statistics for the estimated regression coefficients of the logistic regression model (see logistic regression with IWLS).

The covariance matrix of a logistic logistic regression model with the regression coefficients 𝛼\mathit{\alpha} (with PP coefficients) with a dispersion parameter of 1 of the binomial family is given as:

  • cov̂(𝛼)=(𝐗𝜔)⊤(𝐗𝜔)\hat{\text{cov}}(\mathit{\alpha}) = (\mathbf{X} \mathit{\omega})^\top (\mathbf{X} \mathit{\omega})

where 𝐗\mathbf{X} is the model matrix of the concomitant model (un-standardized) of dimension N×PN \times P and 𝜔=(𝜋*(1−𝜋))12\mathit{\omega} = (\mathit{\pi} * (1 - \mathit{\pi}))^\frac{1}{2}. π\pi is final response, probability returned by the logistic regression model. The diagonal of the covariance matrix contains the variances of the estimated regression coefficients. Thus,

  • sd̂(𝛼)=(diag((𝐗𝜔)⊤(𝐗𝜔)))12\hat{\text{sd}}(\mathit{\alpha}) = \Big( \text{diag}\big((\mathbf{X} \mathit{\omega})^\top (\mathbf{X} \mathit{\omega}) \big) \Big)^\frac{1}{2}.