[BioC] maSigPro - NA values for coefficient estimates - Is a polynomical fit recommended for my data?
jeremy wilson
jeremy.wilson88 at gmail.com
Thu Nov 5 20:05:01 CET 2009
Dear BioConductors,
I have unbalanced number of time points in my experiment. In 3 groups
I have, one group has only observations at one time point (time 0),
the other at two time points (times 0, 5) and the last group with four
time points(0, 5, 12, 31). I am wondering if I can use maSigPro for
this type of data. In this case, what order of polynomial fit is good?
Do you recommend to use degree = 1 to reduce the degree of the model
as the best model is the one with the least order? or is it good to
model with degree = 3? I have 4 time points in one group.
Is it better to do this kind of unbalanced time points data with
LIMMA? There are many comparisons to consider if I use LIMMA. What do
you suggest to use?
I did fit a degree = 3 polynomial model and got NA values. The NA
values appeared as below. There were no NA values for the remaining
coefficients. Are these due to no values for some experimental groups
at some time points? Can I ignore these warnings?
Here are my commands and warnings I got..
design.init<-read.table("maSigPro_PhenotypeInput.txt", sep="\t",
header=TRUE, row.names=1)
> design<-make.design.matrix(design.init, deg=3) # Create a maSigPro design matrix
> design
$dis
VehiclevsNaive BCGvsNaive Time TimexVehicle TimexBCG
Time2 Time2xVehicle Time2xBCG Time3 Time3xVehicle Time3xBCG
Veh Day0 A2.CEL 1 0 0 0 0
0 0 0 0 0 0
Veh Day0 A4.CEL 1 0 0 0 0
0 0 0 0 0 0
Veh Day0 A5.CEL 1 0 0 0 0
0 0 0 0 0 0
Veh Day12 A1.CEL 1 0 12 12 0
144 144 0 1728 1728 0
Veh Day12 A3.CEL 1 0 12 12 0
144 144 0 1728 1728 0
Veh Day0 B1.CEL 1 0 0 0 0
0 0 0 0 0 0
Veh Day0 B3.CEL 1 0 0 0 0
0 0 0 0 0 0
Veh Day0 B4.CEL 1 0 0 0 0
0 0 0 0 0 0
Veh Day12 B2.CEL 1 0 12 12 0
144 144 0 1728 1728 0
Veh Day12 B3.CEL 1 0 12 12 0
144 144 0 1728 1728 0
Veh Day0 C1.CEL 1 0 0 0 0
0 0 0 0 0 0
Veh Day0 C3.CEL 1 0 0 0 0
0 0 0 0 0 0
Veh Day0 C4.CEL 1 0 0 0 0
0 0 0 0 0 0
Veh Day0 C5.CEL 1 0 0 0 0
0 0 0 0 0 0
Veh Day12 C1.CEL 1 0 12 12 0
144 144 0 1728 1728 0
Veh Day12 C3.CEL 1 0 12 12 0
144 144 0 1728 1728 0
Day 0 A1.CEL 0 1 0 0 0
0 0 0 0 0 0
Day 0 A2.CEL 0 1 0 0 0
0 0 0 0 0 0
Day 0 A3.CEL 0 1 0 0 0
0 0 0 0 0 0
Day 5 A1.CEL 0 1 5 0 5
25 0 25 125 0 125
Day 5 A2.CEL 0 1 5 0 5
25 0 25 125 0 125
Day 5 A3.CEL 0 1 5 0 5
25 0 25 125 0 125
Day 12 A1.CEL 0 1 12 0 12
144 0 144 1728 0 1728
Day 12 A2.CEL 0 1 12 0 12
144 0 144 1728 0 1728
Day 12 A3.CEL 0 1 12 0 12
144 0 144 1728 0 1728
Day 31 A1.CEL 0 1 31 0 31
961 0 961 29791 0 29791
Day 31 A2.CEL 0 1 31 0 31
961 0 961 29791 0 29791
Day 31 A3.CEL 0 1 31 0 31
961 0 961 29791 0 29791
Day 0 B1.CEL 0 1 0 0 0
0 0 0 0 0 0
Day 0 B2.CEL 0 1 0 0 0
0 0 0 0 0 0
Day 0 B3.CEL 0 1 0 0 0
0 0 0 0 0 0
Day 5 B1.CEL 0 1 5 0 5
25 0 25 125 0 125
Day 5 B2.CEL 0 1 5 0 5
25 0 25 125 0 125
Day 5 B3.CEL 0 1 5 0 5
25 0 25 125 0 125
Day 12 B1.CEL 0 1 12 0 12
144 0 144 1728 0 1728
Day 12 B2.CEL 0 1 12 0 12
144 0 144 1728 0 1728
Day 12 B3.CEL 0 1 12 0 12
144 0 144 1728 0 1728
Day 31 B1.CEL 0 1 31 0 31
961 0 961 29791 0 29791
Day 31 B2.CEL 0 1 31 0 31
961 0 961 29791 0 29791
Day 31 B3.CEL 0 1 31 0 31
961 0 961 29791 0 29791
Day 0 C1.CEL 0 1 0 0 0
0 0 0 0 0 0
Day 0 C2.CEL 0 1 0 0 0
0 0 0 0 0 0
Day 0 C3.CEL 0 1 0 0 0
0 0 0 0 0 0
Day 0 C4.CEL 0 1 0 0 0
0 0 0 0 0 0
Day 5 C1.CEL 0 1 5 0 5
25 0 25 125 0 125
Day 5 C2.CEL 0 1 5 0 5
25 0 25 125 0 125
Day 5 C3.CEL 0 1 5 0 5
25 0 25 125 0 125
Day 12 C1.CEL 0 1 12 0 12
144 0 144 1728 0 1728
Day 12 C2.CEL 0 1 12 0 12
144 0 144 1728 0 1728
Day 12 C3.CEL 0 1 12 0 12
144 0 144 1728 0 1728
Day 31 C1.CEL 0 1 31 0 31
961 0 961 29791 0 29791
Day 31 C2.CEL 0 1 31 0 31
961 0 961 29791 0 29791
Day 31 C3.CEL 0 1 31 0 31
961 0 961 29791 0 29791
Naive A1.CEL 0 0 0 0 0
0 0 0 0 0 0
Naive A2.CEL 0 0 0 0 0
0 0 0 0 0 0
Naive B1.CEL 0 0 0 0 0
0 0 0 0 0 0
Naive B2.CEL 0 0 0 0 0
0 0 0 0 0 0
Naive C1.CEL 0 0 0 0 0
0 0 0 0 0 0
Naive C2.CEL 0 0 0 0 0
0 0 0 0 0 0
Naive C3.CEL 0 0 0 0 0
0 0 0 0 0 0
$groups.vector
[1] "VehiclevsNaive" "BCGvsNaive" "Naive"
"VehiclevsNaive" "BCGvsNaive" "Naive" "VehiclevsNaive"
"BCGvsNaive"
[9] "Naive" "VehiclevsNaive" "BCGvsNaive"
$edesign
Time Replicate Naive Vehicle BCG
Veh Day0 A2.CEL 0 1 0 1 0
Veh Day0 A4.CEL 0 1 0 1 0
Veh Day0 A5.CEL 0 1 0 1 0
Veh Day12 A1.CEL 12 2 0 1 0
Veh Day12 A3.CEL 12 2 0 1 0
Veh Day0 B1.CEL 0 3 0 1 0
Veh Day0 B3.CEL 0 3 0 1 0
Veh Day0 B4.CEL 0 3 0 1 0
Veh Day12 B2.CEL 12 4 0 1 0
Veh Day12 B3.CEL 12 4 0 1 0
Veh Day0 C1.CEL 0 5 0 1 0
Veh Day0 C3.CEL 0 5 0 1 0
Veh Day0 C4.CEL 0 5 0 1 0
Veh Day0 C5.CEL 0 5 0 1 0
Veh Day12 C1.CEL 12 6 0 1 0
Veh Day12 C3.CEL 12 6 0 1 0
Day 0 A1.CEL 0 7 0 0 1
Day 0 A2.CEL 0 7 0 0 1
Day 0 A3.CEL 0 7 0 0 1
Day 5 A1.CEL 5 8 0 0 1
Day 5 A2.CEL 5 8 0 0 1
Day 5 A3.CEL 5 8 0 0 1
Day 12 A1.CEL 12 9 0 0 1
Day 12 A2.CEL 12 9 0 0 1
Day 12 A3.CEL 12 9 0 0 1
Day 31 A1.CEL 31 10 0 0 1
Day 31 A2.CEL 31 10 0 0 1
Day 31 A3.CEL 31 10 0 0 1
Day 0 B1.CEL 0 11 0 0 1
Day 0 B2.CEL 0 11 0 0 1
Day 0 B3.CEL 0 11 0 0 1
Day 5 B1.CEL 5 12 0 0 1
Day 5 B2.CEL 5 12 0 0 1
Day 5 B3.CEL 5 12 0 0 1
Day 12 B1.CEL 12 13 0 0 1
Day 12 B2.CEL 12 13 0 0 1
Day 12 B3.CEL 12 13 0 0 1
Day 31 B1.CEL 31 14 0 0 1
Day 31 B2.CEL 31 14 0 0 1
Day 31 B3.CEL 31 14 0 0 1
Day 0 C1.CEL 0 15 0 0 1
Day 0 C2.CEL 0 15 0 0 1
Day 0 C3.CEL 0 15 0 0 1
Day 0 C4.CEL 0 15 0 0 1
Day 5 C1.CEL 5 16 0 0 1
Day 5 C2.CEL 5 16 0 0 1
Day 5 C3.CEL 5 16 0 0 1
Day 12 C1.CEL 12 17 0 0 1
Day 12 C2.CEL 12 17 0 0 1
Day 12 C3.CEL 12 17 0 0 1
Day 31 C1.CEL 31 18 0 0 1
Day 31 C2.CEL 31 18 0 0 1
Day 31 C3.CEL 31 18 0 0 1
Naive A1.CEL 0 19 1 0 0
Naive A2.CEL 0 19 1 0 0
Naive B1.CEL 0 20 1 0 0
Naive B2.CEL 0 20 1 0 0
Naive C1.CEL 0 21 1 0 0
Naive C2.CEL 0 21 1 0 0
Naive C3.CEL 0 21 1 0 0
> fit <- p.vector(efiltered.mat, design, Q = 0.05, MT.adjust = "BH",min.obs = 20)
> tstep <- T.fit(fit, step.method = "two.ways.backward", alfa = 0.05)
I got the following warnings after the above step
warnings:
1: In rbind(sol, as.numeric(c(p.value, result$r.squared, ... :
NAs introduced by coercion
2: In rbind(sol, as.numeric(c(p.value, result$r.squared, ... :
NAs introduced by coercion
>sigs <- get.siggenes(tstep, rsq = 0.6, vars = "groups")
> see.genes(sigs$sig.genes$BCGvsNaive, main = "BCGvsNaive", show.fit = T,dis =design$dis, cluster.method="kmeans" ,cluster.data = 1, k = 9)
Warning messages:
1: In rbind(sol, as.numeric(c(p.value, result$r.squared, p.valor))) :
NAs introduced by coercion
2: In rbind(sol, as.numeric(c(p.value, result$r.squared, p.valor))) :
NAs introduced by coercion
betaTimexBCG betaTime2 betaTime2xVehicle
betaTime2xBCG betaTime3
AB002558_at NA 0.0000000000 NA
NA -1.824498e-05
AF001953_g_at NA 0.0114056076 NA
NA -2.616775e-04
AF067795_at NA 0.0011159473 NA
NA -3.539580e-05
D13120_s_at NA 0.0011509775 NA
NA 0.000000e+00
E01534cds_s_at NA -0.0029700508 NA
NA 8.818569e-05
L24896_s_at NA -0.0164766816 NA
NA 3.728959e-04
M20156_at NA -0.0011509235 NA
NA 0.000000e+00
M25638_s_at NA 0.0137020816 NA
NA -3.174337e-04
M58404_at NA -0.0123769559 NA
NA 2.928362e-04
M60921_at NA 0.0000000000 NA
NA -1.465997e-05
M91652complete_seq_at NA -0.0020565596 NA
NA 0.000000e+00
M94918mRNA_f_at NA -0.0256353952 NA
NA 6.500460e-04
M94919mRNA_f_at NA 0.0000000000 NA
NA 3.851441e-05
rc_AA859783_at NA 0.0009030979 NA
NA 0.000000e+00
rc_AA875171_g_at NA -0.0007324529 NA
NA 0.000000e+00
rc_AA891171_s_at NA 0.0010095092 NA
NA 0.000000e+00
rc_AA891727_g_at NA -0.0047867258 NA
NA 1.441130e-04
rc_AA892123_at NA -0.0038206649 NA
NA 1.165693e-04
rc_AA892248_g_at NA 0.0000000000 NA
NA 0.000000e+00
rc_AA893246_at NA 0.0010566522 NA
NA 0.000000e+00
rc_AA924772_at NA -0.0008263709 NA
NA 0.000000e+00
rc_AA944324_at NA -0.0013202929 NA
NA 0.000000e+00
rc_AI007824_g_at NA 0.0000000000 NA
NA 0.000000e+00
rc_AI013194_at NA 0.0000000000 NA
NA 0.000000e+00
rc_AI102044_at NA 0.0000000000 NA
NA 0.000000e+00
rc_AI104544_at NA 0.0019113263 NA
NA 0.000000e+00
rc_AI170268_at NA -0.0012527673 NA
NA 0.000000e+00
rc_AI176662_s_at NA 0.0036907561 NA
NA -1.149829e-04
rc_AI179576_s_at NA -0.0275651514 NA
NA 6.926836e-04
rc_AI235364_at NA -0.0028867693 NA
NA 8.997026e-05
rc_AI237836_at NA 0.0000000000 NA
NA 0.000000e+00
S68135_s_at NA 0.0064503716 NA
NA -1.943265e-04
S74351_s_at NA 0.0066852568 NA
NA -2.153613e-04
U02553cds_s_at NA 0.0000000000 NA
NA -2.903164e-05
U12568_at NA -0.0012220408 NA
NA 0.000000e+00
U19866_at NA -0.0010647292 NA
NA 0.000000e+00
U39875_at NA -0.0007813637 NA
NA 0.000000e+00
U84727_at NA 0.0104752791 NA
NA -2.460059e-04
X03347cds_g_at NA 0.0022984842 NA
NA -7.116160e-05
X54510_g_at NA 0.0000000000 NA
NA 2.572912e-05
X56325mRNA_s_at NA -0.0214234307 NA
NA 5.057676e-04
X58389cds_s_at NA -0.0034907821 NA
NA 1.055627e-04
X63375exon_at NA 0.0132059943 NA
NA -3.082255e-04
X66369_at NA 0.0014723516 NA
NA 0.000000e+00
X94242_at NA 0.0000000000 NA
NA 0.000000e+00
betaTime3xVehicle betaTime3xBCG
AB002558_at NA NA
AF001953_g_at NA NA
AF067795_at NA NA
D13120_s_at NA NA
E01534cds_s_at NA NA
L24896_s_at NA NA
M20156_at NA NA
M25638_s_at NA NA
M58404_at NA NA
M60921_at NA NA
M91652complete_seq_at NA NA
M94918mRNA_f_at NA NA
M94919mRNA_f_at NA NA
rc_AA859783_at NA NA
rc_AA875171_g_at NA NA
rc_AA891171_s_at NA NA
rc_AA891727_g_at NA NA
rc_AA892123_at NA NA
rc_AA892248_g_at NA NA
rc_AA893246_at NA NA
rc_AA924772_at NA NA
rc_AA944324_at NA NA
rc_AI007824_g_at NA NA
rc_AI013194_at NA NA
rc_AI102044_at NA NA
rc_AI104544_at NA NA
rc_AI170268_at NA NA
rc_AI176662_s_at NA NA
rc_AI179576_s_at NA NA
rc_AI235364_at NA NA
rc_AI237836_at NA NA
S68135_s_at NA NA
S74351_s_at NA NA
U02553cds_s_at NA NA
U12568_at NA NA
U19866_at NA NA
U39875_at NA NA
U84727_at NA NA
X03347cds_g_at NA NA
X54510_g_at NA NA
X56325mRNA_s_at NA NA
X58389cds_s_at NA NA
X63375exon_at NA NA
X66369_at NA NA
X94242_at NA NA
The following coefficients had no "NA" values.
tstep$coefficients$beta0
tstep$coefficients$betaVehiclevsNaive
tstep$coefficients$betaBCGvsNaive
tstep$coefficients$betaTime tstep$coefficients$betaTimexVehicle
> sessionInfo()
R version 2.10.0 (2009-10-26)
i386-pc-mingw32
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] tcltk stats graphics grDevices utils datasets
methods base
other attached packages:
[1] maSigPro_1.17.5 DynDoc_1.23.1
widgetTools_1.23.1 genefilter_1.26.4
[5] arrayQualityMetrics_2.4.1 affyPLM_1.22.0
preprocessCore_1.7.9 gcrma_2.17.4
[9] rgu34acdf_2.5.0 affy_1.23.12 Biobase_2.5.8
loaded via a namespace (and not attached):
[1] affyio_1.13.5 annotate_1.23.4 AnnotationDbi_1.7.20
beadarray_1.13.9
[5] Biostrings_2.13.54 DBI_0.2-4 grid_2.10.0
hwriter_1.1
[9] IRanges_1.3.99 KernSmooth_2.23-3 lattice_0.17-26
latticeExtra_0.6-3
[13] limma_3.0.3 marray_1.23.0 Mfuzz_2.3.4
RColorBrewer_1.0-2
[17] RSQLite_0.7-3 simpleaffy_2.21.3 splines_2.10.0
stats4_2.10.0
[21] survival_2.35-7 tkWidgets_1.23.2 tools_2.10.0
vsn_3.13.7
[25] xtable_1.5-5
I would really appreciate for any suggestions. Thanks in advance..
More information about the Bioconductor
mailing list