summary function does not work with Westfall correction in multcomp with 27 comparisons

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

summary function does not work with Westfall correction in multcomp with 27 comparisons

Friedman, Richard A.
Dear List.

I ran multcomp with 27 comaprisons. The glht command returned an mcp object,
but the summary command with the Westfall correction ddi not give a summary.
When I ran the same dataset with 4 comparisons I got p-values. When I sued a summary with
univariate or Bonferroni’s method with all 27 comarisons I got p-values. But all 27 did not
work for Wesrfall. Please advise.
Here is a record of my session with 27 comparisons and Westfall:

R version 3.5.1 (2018-07-02) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin15.6.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[R.app GUI 1.70 (7543) x86_64-apple-darwin15.6.0]

[History restored from /Users/friedman/.Rapp.history]

objc[31790]: Class FIFinderSyncExtensionHost is implemented in both /System/Library/PrivateFrameworks/FinderKit.framework/Versions/A/FinderKit (0x7fff9c2b7c90) and /System/Library/PrivateFrameworks/FileProvider.framework/OverrideBundles/FinderSyncCollaborationFileProviderOverride.bundle/Contents/MacOS/FinderSyncCollaborationFileProviderOverride (0x119bf8cd8). One of the two will be used. Which one is undefined.
> library(multcomp)
Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data
Loading required package: MASS

Attaching package: ‘TH.data’

The following object is masked from ‘package:MASS’:

    geyser

> tumor<-read.table("a8_1wayall_input.txt",sep="\t",header=T)
> class(tumor)
[1] "data.frame"
> dim(tumor)
[1] 309   2
> tumor[1,]
     condition  log2vol
1 aa.vector.4w 7.297375
>
> model<-lm(log2vol~condition,data=tumor)
> summary(model)

Call:
lm(formula = log2vol ~ condition, data = tumor)

Residuals:
   Min     1Q Median     3Q    Max
-7.997 -2.414  0.291  2.164  8.059

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)
(Intercept)                   4.58776    0.60457   7.588  4.3e-13 ***
conditionab.vector.6w         1.39818    0.85499   1.635 0.103054
conditionac.vector.8w         2.89085    0.85499   3.381 0.000820 ***
conditionba.dnmt1kd.4w       -2.67491    0.84733  -3.157 0.001760 **
conditionbb.dnmt1kd.6w       -1.43390    0.84733  -1.692 0.091654 .
conditionbc.dnmt1kd.8w       -0.47188    0.84733  -0.557 0.578020
conditionca.dnmt3bkd.4w      -3.15325    0.89139  -3.537 0.000469 ***
conditioncb.dnmt3bkd.6w      -2.17334    0.89139  -2.438 0.015355 *
conditioncc.dnmt3bkd.8w      -1.50187    0.89139  -1.685 0.093078 .
conditionda.dnmt1kdhrad9.4w   1.08153    1.08991   0.992 0.321863
conditiondb.dnmt1kdhrad9.6w   2.38383    1.08991   2.187 0.029517 *
conditiondc.dnmt1kdhrad9.8w   3.53990    1.08991   3.248 0.001297 **
conditionea.dnmt3bkdhrad9.4w  0.02795    1.06049   0.026 0.978991
conditioneb.dnmt3bkdhrad9.6w  1.34037    1.06049   1.264 0.207263
conditionec.dnmt3bkdhrad9.8w  3.40955    1.06049   3.215 0.001449 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.141 on 294 degrees of freedom
Multiple R-squared:  0.3155, Adjusted R-squared:  0.2829
F-statistic: 9.678 on 14 and 294 DF,  p-value: < 2.2e-16

> model.mc<-glht(model,linfct=mcp(condition=c("ab.vector.6w-aa.vector.4w=0",
+                        "ac.vector.8w-aa.vector.4w=0",
+                        "ac.vector.8w-ab.vector.6w=0",
+       "bb.dnmt1kd.6w-ba.dnmt1kd.4w=0",
+                        "bc.dnmt1kd.8w-ba.dnmt1kd.4w=0",
+                        "bc.dnmt1kd.8w-bb.dnmt1kd.6w=0",
+         "cb.dnmt3bkd.6w-ca.dnmt3bkd.4w=0",
+                        "cc.dnmt3bkd.8w-ca.dnmt3bkd.4w=0",
+                        "cc.dnmt3bkd.8w-cb.dnmt3bkd.6w=0",
+         "db.dnmt1kdhrad9.6w-da.dnmt1kdhrad9.4w=0",
+                        "dc.dnmt1kdhrad9.8w-da.dnmt1kdhrad9.4w=0",
+                        "dc.dnmt1kdhrad9.8w-db.dnmt1kdhrad9.6w=0",
+         "eb.dnmt3bkdhrad9.6w-ea.dnmt3bkdhrad9.4w=0",
+                        "ec.dnmt3bkdhrad9.8w-ea.dnmt3bkdhrad9.4w=0",
+                        "ec.dnmt3bkdhrad9.8w-eb.dnmt3bkdhrad9.6w=0",
+       "ba.dnmt1kd.4w-aa.vector.4w=0",
+                        "ca.dnmt3bkd.4w-aa.vector.4w=0",
+                        "da.dnmt1kdhrad9.4w-ba.dnmt1kd.4w=0",
+                        "ea.dnmt3bkdhrad9.4w-ca.dnmt3bkd.4w=0",
+       "bb.dnmt1kd.6w-ab.vector.6w=0",
+                        "cb.dnmt3bkd.6w-ab.vector.6w=0",
+                        "db.dnmt1kdhrad9.6w-bb.dnmt1kd.6w=0",
+                        "eb.dnmt3bkdhrad9.6w-cb.dnmt3bkd.6w=0",
+       "bc.dnmt1kd.8w-ac.vector.8w=0",
+                        "cc.dnmt3bkd.8w-ac.vector.8w=0",
+                        "dc.dnmt1kdhrad9.8w-bc.dnmt1kd.8w=0",
+                        "ec.dnmt3bkdhrad9.8w-cc.dnmt3bkd.8w=0")))
>
> summary(model.mc,test=adjusted(type="Westfall"))

######################
THE PROGRAM FROZE AT THIS POINT AND DID NOT RETURN ANYTHING
HERE IS MY SESSION INFO

######################

> sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] multcomp_1.4-8  TH.data_1.0-9   MASS_7.3-50     survival_2.42-6 mvtnorm_1.0-8

loaded via a namespace (and not attached):
[1] zoo_1.8-3        compiler_3.5.1   Matrix_1.2-14    sandwich_2.4-0   codetools_0.2-15 splines_3.5.1    grid_3.5.1       lattice_0.20-35

#############################

Why would the summary function not work with 27 comparisons with Westfall but work with 4 comparisons?
Why would the summary function not work with 27 comparisons with Westfall but work with the same  comparisons with Bonferroni?
Is the problem intrinsic to
A. Westfall’s , method with a large number of comarisons.
B, The impkementation of Westfall’s method in multcomp?
C. Machine requirements for Westfall’s method with so many comparisons?
D. My coding
E. Other.

I have used Westsall;’s method throught the paper which I am now working and would prefer to use it for this
problem for consistency.
I would appreciate any advice.


Thanks and best wishes,
Rich
Richard A. Friedman, PhD
Associate Research Scientist,
Biomedical Informatics Shared Resource
Herbert Irving Comprehensive Cancer Center (HICCC)
Lecturer,
Department of Biomedical Informatics (DBMI)
Room 825
Irving Cancer Research Center
Columbia University Herbert and Florence Irving Medical Center
1130 St. Nicholas Ave
New York, NY 10032
(212)851-4765 (voice)
[hidden email]<mailto:[hidden email]>

http://www.columbia.edu/~raf4/index.html


In memoriam, Steve Ditko


        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: summary function does not work with Westfall correction in multcomp with 27 comparisons

Gerrit Eichner
Dear Rich,

w/o the original data we can actually only guess, but I think you
may haven't been patient enough to let summary.glht finish its
job. Have you tried to increase the number of contrasts, i.e.
comparisons, step by step to see how the computational burden and
hence the required computing time increases?

Internally, a lot effort goes into the computation of probabilities
and/or quantiles of multivariate normal or t distributions, and in
your setting they *are* high-dimensional. For more reliable and
elaborate details you may have to consult the cited references in
?summary.glht, or wait/hope for a more knowledgeable list member to
"jump in".

  Hth --  Gerrit

---------------------------------------------------------------------
Dr. Gerrit Eichner                   Mathematical Institute, Room 212
[hidden email]   Justus-Liebig-University Giessen
Tel: +49-(0)641-99-32104          Arndtstr. 2, 35392 Giessen, Germany
http://www.uni-giessen.de/eichner
---------------------------------------------------------------------

Am 08.11.2018 um 18:10 schrieb Friedman, Richard A.:

> Dear List.
>
> I ran multcomp with 27 comaprisons. The glht command returned an mcp object,
> but the summary command with the Westfall correction ddi not give a summary.
> When I ran the same dataset with 4 comparisons I got p-values. When I sued a summary with
> univariate or Bonferroni’s method with all 27 comarisons I got p-values. But all 27 did not
> work for Wesrfall. Please advise.
> Here is a record of my session with 27 comparisons and Westfall:
>
> R version 3.5.1 (2018-07-02) -- "Feather Spray"
> Copyright (C) 2018 The R Foundation for Statistical Computing
> Platform: x86_64-apple-darwin15.6.0 (64-bit)
>
> R is free software and comes with ABSOLUTELY NO WARRANTY.
> You are welcome to redistribute it under certain conditions.
> Type 'license()' or 'licence()' for distribution details.
>
>    Natural language support but running in an English locale
>
> R is a collaborative project with many contributors.
> Type 'contributors()' for more information and
> 'citation()' on how to cite R or R packages in publications.
>
> Type 'demo()' for some demos, 'help()' for on-line help, or
> 'help.start()' for an HTML browser interface to help.
> Type 'q()' to quit R.
>
> [R.app GUI 1.70 (7543) x86_64-apple-darwin15.6.0]
>
> [History restored from /Users/friedman/.Rapp.history]
>
> objc[31790]: Class FIFinderSyncExtensionHost is implemented in both /System/Library/PrivateFrameworks/FinderKit.framework/Versions/A/FinderKit (0x7fff9c2b7c90) and /System/Library/PrivateFrameworks/FileProvider.framework/OverrideBundles/FinderSyncCollaborationFileProviderOverride.bundle/Contents/MacOS/FinderSyncCollaborationFileProviderOverride (0x119bf8cd8). One of the two will be used. Which one is undefined.
>> library(multcomp)
> Loading required package: mvtnorm
> Loading required package: survival
> Loading required package: TH.data
> Loading required package: MASS
>
> Attaching package: ‘TH.data’
>
> The following object is masked from ‘package:MASS’:
>
>      geyser
>
>> tumor<-read.table("a8_1wayall_input.txt",sep="\t",header=T)
>> class(tumor)
> [1] "data.frame"
>> dim(tumor)
> [1] 309   2
>> tumor[1,]
>       condition  log2vol
> 1 aa.vector.4w 7.297375
>>
>> model<-lm(log2vol~condition,data=tumor)
>> summary(model)
>
> Call:
> lm(formula = log2vol ~ condition, data = tumor)
>
> Residuals:
>     Min     1Q Median     3Q    Max
> -7.997 -2.414  0.291  2.164  8.059
>
> Coefficients:
>                               Estimate Std. Error t value Pr(>|t|)
> (Intercept)                   4.58776    0.60457   7.588  4.3e-13 ***
> conditionab.vector.6w         1.39818    0.85499   1.635 0.103054
> conditionac.vector.8w         2.89085    0.85499   3.381 0.000820 ***
> conditionba.dnmt1kd.4w       -2.67491    0.84733  -3.157 0.001760 **
> conditionbb.dnmt1kd.6w       -1.43390    0.84733  -1.692 0.091654 .
> conditionbc.dnmt1kd.8w       -0.47188    0.84733  -0.557 0.578020
> conditionca.dnmt3bkd.4w      -3.15325    0.89139  -3.537 0.000469 ***
> conditioncb.dnmt3bkd.6w      -2.17334    0.89139  -2.438 0.015355 *
> conditioncc.dnmt3bkd.8w      -1.50187    0.89139  -1.685 0.093078 .
> conditionda.dnmt1kdhrad9.4w   1.08153    1.08991   0.992 0.321863
> conditiondb.dnmt1kdhrad9.6w   2.38383    1.08991   2.187 0.029517 *
> conditiondc.dnmt1kdhrad9.8w   3.53990    1.08991   3.248 0.001297 **
> conditionea.dnmt3bkdhrad9.4w  0.02795    1.06049   0.026 0.978991
> conditioneb.dnmt3bkdhrad9.6w  1.34037    1.06049   1.264 0.207263
> conditionec.dnmt3bkdhrad9.8w  3.40955    1.06049   3.215 0.001449 **
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Residual standard error: 3.141 on 294 degrees of freedom
> Multiple R-squared:  0.3155, Adjusted R-squared:  0.2829
> F-statistic: 9.678 on 14 and 294 DF,  p-value: < 2.2e-16
>
>> model.mc<-glht(model,linfct=mcp(condition=c("ab.vector.6w-aa.vector.4w=0",
> +                        "ac.vector.8w-aa.vector.4w=0",
> +                        "ac.vector.8w-ab.vector.6w=0",
> +       "bb.dnmt1kd.6w-ba.dnmt1kd.4w=0",
> +                        "bc.dnmt1kd.8w-ba.dnmt1kd.4w=0",
> +                        "bc.dnmt1kd.8w-bb.dnmt1kd.6w=0",
> +         "cb.dnmt3bkd.6w-ca.dnmt3bkd.4w=0",
> +                        "cc.dnmt3bkd.8w-ca.dnmt3bkd.4w=0",
> +                        "cc.dnmt3bkd.8w-cb.dnmt3bkd.6w=0",
> +         "db.dnmt1kdhrad9.6w-da.dnmt1kdhrad9.4w=0",
> +                        "dc.dnmt1kdhrad9.8w-da.dnmt1kdhrad9.4w=0",
> +                        "dc.dnmt1kdhrad9.8w-db.dnmt1kdhrad9.6w=0",
> +         "eb.dnmt3bkdhrad9.6w-ea.dnmt3bkdhrad9.4w=0",
> +                        "ec.dnmt3bkdhrad9.8w-ea.dnmt3bkdhrad9.4w=0",
> +                        "ec.dnmt3bkdhrad9.8w-eb.dnmt3bkdhrad9.6w=0",
> +       "ba.dnmt1kd.4w-aa.vector.4w=0",
> +                        "ca.dnmt3bkd.4w-aa.vector.4w=0",
> +                        "da.dnmt1kdhrad9.4w-ba.dnmt1kd.4w=0",
> +                        "ea.dnmt3bkdhrad9.4w-ca.dnmt3bkd.4w=0",
> +       "bb.dnmt1kd.6w-ab.vector.6w=0",
> +                        "cb.dnmt3bkd.6w-ab.vector.6w=0",
> +                        "db.dnmt1kdhrad9.6w-bb.dnmt1kd.6w=0",
> +                        "eb.dnmt3bkdhrad9.6w-cb.dnmt3bkd.6w=0",
> +       "bc.dnmt1kd.8w-ac.vector.8w=0",
> +                        "cc.dnmt3bkd.8w-ac.vector.8w=0",
> +                        "dc.dnmt1kdhrad9.8w-bc.dnmt1kd.8w=0",
> +                        "ec.dnmt3bkdhrad9.8w-cc.dnmt3bkd.8w=0")))
>>
>> summary(model.mc,test=adjusted(type="Westfall"))
>
> ######################
> THE PROGRAM FROZE AT THIS POINT AND DID NOT RETURN ANYTHING
> HERE IS MY SESSION INFO
>
> ######################
>
>> sessionInfo()
> R version 3.5.1 (2018-07-02)
> Platform: x86_64-apple-darwin15.6.0 (64-bit)
> Running under: macOS High Sierra 10.13.6
>
> Matrix products: default
> BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
> LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] multcomp_1.4-8  TH.data_1.0-9   MASS_7.3-50     survival_2.42-6 mvtnorm_1.0-8
>
> loaded via a namespace (and not attached):
> [1] zoo_1.8-3        compiler_3.5.1   Matrix_1.2-14    sandwich_2.4-0   codetools_0.2-15 splines_3.5.1    grid_3.5.1       lattice_0.20-35
>
> #############################
>
> Why would the summary function not work with 27 comparisons with Westfall but work with 4 comparisons?
> Why would the summary function not work with 27 comparisons with Westfall but work with the same  comparisons with Bonferroni?
> Is the problem intrinsic to
> A. Westfall’s , method with a large number of comarisons.
> B, The impkementation of Westfall’s method in multcomp?
> C. Machine requirements for Westfall’s method with so many comparisons?
> D. My coding
> E. Other.
>
> I have used Westsall;’s method throught the paper which I am now working and would prefer to use it for this
> problem for consistency.
> I would appreciate any advice.
>
>
> Thanks and best wishes,
> Rich
> Richard A. Friedman, PhD
> Associate Research Scientist,
> Biomedical Informatics Shared Resource
> Herbert Irving Comprehensive Cancer Center (HICCC)
> Lecturer,
> Department of Biomedical Informatics (DBMI)
> Room 825
> Irving Cancer Research Center
> Columbia University Herbert and Florence Irving Medical Center
> 1130 St. Nicholas Ave
> New York, NY 10032
> (212)851-4765 (voice)
> [hidden email]<mailto:[hidden email]>
>
> http://www.columbia.edu/~raf4/index.html
>
>
> In memoriam, Steve Ditko
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.