Introduction Simultaneous Inference in General Parametric Models Torsten Hothorn Institut f¨ ur Statistik (in collaboration with Frank Bretz, Novartis, and Peter Westfall, Texas Tech) 1 WU Wien, 2009-01-23 Introduction Introduction 2 3 WU Wien, 2009-01-23 WU Wien, 2009-01-23
Introduction Model M (( Z 1 , . . . , Z n ) , θ, η ) is a (semiparametric) model with • n observations ( Z 1 , . . . , Z n ) • elemental parameters θ ∈ R p and • other (random or nuisance) parameters η . We are interested in linear functions ϑ := K θ defined by a constant matrix K ∈ R k,p . 4 5 WU Wien, 2009-01-23 WU Wien, 2009-01-23 Estimation Distribution of ϑ θ n ∈ R p is an estimate of θ and S n ∈ R p,p is an estimate of cov(ˆ ˆ θ n ) with By Theorem 3.3.A in Serfling (1980), the linear function ˆ ϑ n = K ˆ θ n , i.e., P → Σ ∈ R p,p a n S n − an estimate of our parameters of interest, also follows an approximate multivariate normal distribution for some positive, nondecreasing sequence a n . a ∼ N k ( ϑ, S ⋆ ϑ n = K ˆ ˆ n ) θ n A multivariate central limit theorem is assumed: n := K S n K ⊤ for any fixed matrix K ∈ R k,p with covariance matrix S ⋆ a 1 / 2 d (ˆ θ n − θ ) − → N p (0 , Σ) . n Therefore, we simply assume a We write ˆ a P → Σ ⋆ := K Σ K ⊤ ∈ R k,k ∼ N p ( θ, S n ). ∼ N k ( ϑ, S ⋆ n ) with a n S ⋆ θ n ˆ ϑ n − n These assumptions are fulfilled for most of the models commonly in use. 6 7 WU Wien, 2009-01-23 WU Wien, 2009-01-23
A Statistic and its Distribution General Linear Hypothesis Consider the multivariate statistic Consider the null hypothesis T n := D − 1 / 2 (ˆ ϑ n − ϑ ) n H 0 : ϑ := K θ = m . where D n = diag(S ⋆ n ) is the diagonal matrix given by the diagonal elements of S ⋆ n . Classically, F - or χ 2 -statistics are used to test H 0 . However, a rejection of H 0 does not give further indication about the nature of the significant By Slutzky’s Theorem, this statistic is again asymptotically normally result. Therefore, one is often interested in the individual null hypotheses distributed H j 0 : ϑ j = m j . a ∼ N k (0 , R n ) T n where Testing the hypotheses set { H 1 0 , . . . , H k 0 } simultaneously thus requires the R n = D − 1 / 2 n D − 1 / 2 S ⋆ ∈ R k,k individual assessments while maintaining the familywise error rate. n n is the correlation matrix of the k -dimensional statistic T n . 8 9 WU Wien, 2009-01-23 WU Wien, 2009-01-23 A Maximum-Type Statistic Null-Distribution and a Global Test t t � � P (max( | T n | ) ≤ t ) ∼ · · · ϕ k ( x 1 , . . . , x k ; R ) dx 1 · · · dx k =: g ( R , t ) An alternative test statistic for testing H 0 is = − t − t max( | T n | ) where ϕ k is the k -dimensional normal density function. Can we approximate it’s distribution under H 0 efficiently? R is not known but g ( R , t ) is a continuous function of R and converges P as R n − → R . The integral can be approximated by quasi-randomized We have to find a good approximation of P (max( | T n | ) ≤ t ) for some Monte-Carlo methods (Genz, 1992, Genz and Bretz, 1999). t ∈ R + . The resulting global p -value for H 0 is then p global = 1 − g ( R n , max | t | ) when T = t has been observed. 10 11 WU Wien, 2009-01-23 WU Wien, 2009-01-23
Simultaneous Inference Simultaneous Confidence Intervals But what about the partial hypotheses H 1 0 , . . . , H k 0 ? It’s simple! A simultaneous (1 − 2 α ) × 100% confidence interval for ϑ is given by The multiplicity adjusted p -value for the j th individual two-sided ϑ n ± q α diag( D n ) 1 / 2 ˆ hypothesis where q α is the (approximate) 1 − α quantile of the distribution of H j 0 : ϑ j = m j , j = 1 , . . . , k, max( | T n | ). is given by p j = 1 − g ( R n , | t j | ) , where t = ( t 1 , . . . , t k ) denote the observed test statistics (single-step procedure). Reject each H j 0 at familywise error rate α when p j ≤ α . 12 13 WU Wien, 2009-01-23 WU Wien, 2009-01-23 Examples: Linear Regression Predicting Body Fat Z i = ( Y i , X i ) , i = 1 , . . . , n , with response Y i and exploratory variables X i = ( X i 1 , . . . , X iq ) Garcia et al. (2005) describe a linear model for total body fat prediction. Model: q � Y i = β 0 + β j X ij + σε i , Aim: Based on p = 9 simple measurements (circumferences of elbow, j =1 knee etc) we want to estimate a simple (!) formula to predict the total with elemental parameters θ = ( β 0 , β 1 , . . . , β q ) estimated via body fat obtained for n = 71 healthy German women by means of Dual � − 1 X ⊤ Y ∼ N q +1 � � − 1 � Energy X-Ray Absorptiometry. � X ⊤ X θ, σ 2 � X ⊤ X ˆ θ n = . Problem: Variable selection! Now � − 1 K ⊤ ) θ n ∼ N k ( K θ, σ 2 K � X ⊤ X ϑ n = K ˆ ˆ and T n = D − 1 / 2 ˆ ϑ n ∼ t q +1 ( n − q, R ) exact inference possible! n 14 15 WU Wien, 2009-01-23 WU Wien, 2009-01-23
Linear Model Fit Parameters of Interest R> data("bodyfat", package = "mboost") R> lmod <- lm(DEXfat ~ ., data = bodyfat) R> library("multcomp") R> summary(lmod) R> K <- cbind(0, diag(length(coef(lmod)) - 1)) R> rownames(K) <- names(coef(lmod))[-1] Estimate Std. Error t value Pr(>|t|) R> lmod_glht <- glht(lmod, linfct = K) (Intercept) -69.028276 7.516860 -9.1831 4.184e-13 *** R> K age 0.019962 0.032213 0.6197 0.537767 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] waistcirc 0.210487 0.067145 3.1348 0.002644 ** age 0 1 0 0 0 0 0 0 0 0 hipcirc 0.343513 0.080373 4.2740 6.852e-05 *** waistcirc 0 0 1 0 0 0 0 0 0 0 elbowbreadth -0.412369 1.022907 -0.4031 0.688259 hipcirc 0 0 0 1 0 0 0 0 0 0 kneebreadth 1.757984 0.724952 2.4250 0.018286 * elbowbreadth 0 0 0 0 1 0 0 0 0 0 anthro3a 5.742295 5.207524 1.1027 0.274492 kneebreadth 0 0 0 0 0 1 0 0 0 0 anthro3b 9.866431 5.657864 1.7438 0.086224 . anthro3a 0 0 0 0 0 0 1 0 0 0 anthro3c 0.387430 2.087463 0.1856 0.853376 anthro3b 0 0 0 0 0 0 0 1 0 0 anthro4 -6.574395 6.489177 -1.0131 0.314999 anthro3c 0 0 0 0 0 0 0 0 1 0 --- anthro4 0 0 0 0 0 0 0 0 0 1 Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Multiple R-squared: 0.923, Adjusted R-squared: 0.912 F-statistic: 81.3 on 9 and 61 DF, p-value: <2e-16 16 17 WU Wien, 2009-01-23 WU Wien, 2009-01-23 F -Test Maximum Test R> summary(lmod_glht) R> summary(lmod_glht, test = Ftest()) Simultaneous Tests for General Linear Hypotheses General Linear Hypotheses Fit: lm(formula = DEXfat ~ ., data = bodyfat) Linear Hypotheses: Estimate Linear Hypotheses: age == 0 0.01996 Estimate Std. Error t value Pr(>|t|) waistcirc == 0 0.21049 age == 0 0.01996 0.03221 0.620 0.9959 hipcirc == 0 0.34351 waistcirc == 0 0.21049 0.06714 3.135 0.0213 * elbowbreadth == 0 -0.41237 hipcirc == 0 0.34351 0.08037 4.274 <0.001 *** kneebreadth == 0 1.75798 elbowbreadth == 0 -0.41237 1.02291 -0.403 0.9998 anthro3a == 0 5.74230 kneebreadth == 0 1.75798 0.72495 2.425 0.1316 anthro3b == 0 9.86643 anthro3a == 0 5.74230 5.20752 1.103 0.8948 anthro3c == 0 0.38743 anthro3b == 0 9.86643 5.65786 1.744 0.4783 anthro4 == 0 -6.57439 anthro3c == 0 0.38743 2.08746 0.186 1.0000 anthro4 == 0 -6.57439 6.48918 -1.013 0.9295 Global Test: --- F DF1 DF2 Pr(>F) Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 1 81.35 9 61 1.387e-30 (Adjusted p values reported -- single-step method) 18 19 WU Wien, 2009-01-23 WU Wien, 2009-01-23
ANOVA Genetic Components of Alcoholism Model: n = 24 n = 58 n = 15 Y ij = µ + γ j + ε ij 6 ● ● Overparameterized, usually the elemental parameters are θ = ( µ, γ 2 − γ 1 , γ 3 − γ 1 , . . . , γ q − γ 1 ). 4 Expression Level Dunnett many-to-one comparisons : 2 = (0 , diag( q )) K Dunnett 0 ϑ Dunnett = K Dunnett θ = ( γ 2 − γ 1 , γ 3 − γ 1 , . . . , γ q − γ 1 ) −2 ● Tukey all-pair comparisons : ● ● ● 0 1 0 short intermediate long K Tukey = 0 0 1 NACP−REP1 Allele Length 0 1 − 1 ϑ Tukey = K Tukey θ = ( γ 2 − γ 1 , γ 3 − γ 1 , γ 2 − γ 3 ) 20 21 WU Wien, 2009-01-23 WU Wien, 2009-01-23 Genetic Components of Alcoholism Genetic Components of Alcoholism R> data("alpha", package = "coin") R> amod_glht_sw <- glht(amod, linfct = mcp(alength = "Tukey"), R> amod <- aov(elevel ~ alength, data = alpha) + vcov = sandwich) R> confint(glht(amod, linfct = mcp(alength = "Tukey"))) R> confint(amod_glht_sw) Simultaneous Confidence Intervals Simultaneous Confidence Intervals Multiple Comparisons of Means: Tukey Contrasts Multiple Comparisons of Means: Tukey Contrasts Fit: aov(formula = elevel ~ alength, data = alpha) Fit: aov(formula = elevel ~ alength, data = alpha) Estimated Quantile = 2.3714 Estimated Quantile = 2.3718 95% family-wise confidence level 95% family-wise confidence level Linear Hypotheses: Linear Hypotheses: Estimate lwr upr Estimate lwr upr intermediate - short == 0 0.43415 -0.47561 1.34391 intermediate - short == 0 0.4341523 -0.5713432 1.4396478 long - short == 0 1.18875 -0.04498 2.42248 long - short == 0 1.1887500 0.1376593 2.2398407 long - intermediate == 0 0.75460 -0.33118 1.84038 long - intermediate == 0 0.7545977 -0.0005049 1.5097003 22 23 WU Wien, 2009-01-23 WU Wien, 2009-01-23
Recommend
More recommend