2017-11-26 33 views
-2

Bonferroni補正用のRStudioとPearson Correlation Matrix用の生のP値のコーディングを続行するにはどうすればいいですか?私は学生でRの新人です。データの平均、SD、およびnの表をどのように取得するかについても迷っています。私がピアソン相関行列を計算したとき、私はちょうどr値を得て、生の確率値も得ませんでした。私はどのようにRStudioでそれを得るためにコード化するのか分かりません。私はBonferroni Correctionを計算しようとしましたが、リストオブジェクトはdouble型に強制することはできないというエラーメッセージを受け取りました。これがなくなるようにコードを修正するにはどうすればよいですか?私もデータの平均、SD、nの表を作成しようとしました。R:Bonferroni補正のエラーコードを修正する方法

次のように私のデータは次のとおりです。私はdplyr、tidyr、multcomp、multcompview 私はエクセルCSVからデータを読み込み(カンマ区切り)ファイルと、これは、データを作成しインストールしている

Tree Height DBA Leaf Diameter 
45.3 14.9 0.76 
75.2 26.6 1.06 
70.1 22.9 1.19 
95  31.8 1.59 
107.8 35.5 0.43 
93  26.2 1.49 
91.5 29  1.19 
78.5 29.2 1.1 
85.2 30.3 1.24 
50  16.8 0.67 
47.1 12.8 0.98 
73.2 28.4 1.2 

パッケージ>をdataHW8_1 12obs。 3つの変数

summary(dataHW8_1) 

の私は、データ プロット、xlab = "樹高(dataHW8_1 $ Tree_Height、dataHW8_1 $ DBA、メイン= "バストの高さ(DBA)のトランク直径対散布樹高" の散布図を作成しました(cm)」、「ylab = "DBA(cm)") plot(dataHW8_1 $ Tree_Height、dataHW8_1 $ Leaf_Diameter、main = "散布図ツリー高さvsリーフ直径"、xlab = "ツリー高さ(cm)"、ylab = "リーフ(DBA)と葉の直径での散布図の幹直径、xlab = "DBA(cm)"、ylab = "葉の直径(cm)") plot(dataHW8_1 $ DBA、dataHW8_1 $ Leaf_Diameter、main = )」)

私はデータが線形ではないことに気付きましたrは私は

私は、私はその後、散布

のマトリックスプロットを作成

plot(dataHW8_1log$Tree_Height,dataHW8_1log$DBA,main="Scatterplot of 
Transformed (log)Tree Height Vs Trunk Diameter at Breast Height 
(DBA)",xlab="Tree Height (cm)",ylab="DBA (cm)") 
plot(dataHW8_1log$Tree_Height,dataHW8_1log$Leaf_Diameter,main="Scatterplot 
of Transformed (log)Tree Height Vs Leaf Diameter",xlab="Tree Height 
(cm)",ylab="Leaf Diameter (cm)") 
plot(dataHW8_1log$DBA,dataHW8_1log$Leaf_Diameter,main="Scatterplot of 
Transformed (log) Trunk Diameter at Breast Height (DBA) Vs Leaf 
Diameter",xlab="DBA (cm)",ylab="Leaf Diameter (cm)") 

変換されたデータを使用して私の散布を再作成ログ()fucntion dataHW8_1log =ログ(dataHW8_1)を使用して形質転換されました

pairs(dataHW8_1log) 

次に、Pearson法を使用して相関係数を計算しました これはP値の未補正行列を与えません------どのようにしますか?あなたはそれを行うのですか

cor(dataHW8_1log,method="pearson") 

私はその後、ボンフェローニ補正を計算

データの生確率の行列(未補正P値)を取得するために何をすべきかに引っかかっています-----どのように?

これを行う
Data$Bonferroni = 
p.adjust(dataHW8_1log, 
method = "bonferroni") 

は私にfollwingエラーました:

Error in p.adjust(dataHW8_1log, method = "bonferroni") : 
(list) object cannot be coerced to type 'double' 

を私はlapply使用して修正しようとしたが、それは私のpromblem

を固定していなかった私は、平均値の表を作ってみました、 SD、n、しかし、私は次のコードを作成することしかできず、そこからどこに行くべきかに固執していました.-あなたはどうやってそれをしますか?

(,data = dataHW8_1log, 
     FUN = function(x) c(Mean = mean(x, na.rm = T), 
          n = length(x), 
          sd = sd(x, na.rm = T)) 

私はオンライン次の例を試してみましたが、それらのどれも正しくコードにボンフェローニ補正を得ることで私を助けていません。誰かが私が間違ったことを説明し、マトリックス/テーブルを作る方法を教えてもらえると大変感謝しています。

答えて

0

ここでは、50行×10列のサンプルデータフレームを使用する例を示します。

対散布図を示すことができます。

# Pairwise scatterplot 
pairs(df); 

enter image description here

現在、単一の比較のためのp値を取得するために cor.testを使用することができます。便利な関数 cor.test.pを使用して、すべてのペアワイズ比較でこれを行います。信用供与期限が到来する場合には、 cor.test.p関数は this SO postから取得され、訂正されていないp値の行列を返す間に引数として dataframeをとります。

# cor.test on dataframes 
# From: https://stackoverflow.com/questions/13112238/a-matrix-version-of-cor-test 
cor.test.p <- function(x) { 
    FUN <- function(x, y) cor.test(x, y)[["p.value"]]; 
    z <- outer(
     colnames(x), 
     colnames(x), 
     Vectorize(function(i,j) FUN(x[,i], x[,j]))); 
    dimnames(z) <- list(colnames(x), colnames(x)); 
    return(z); 
} 

# Uncorrected p-values from pairwise correlation tests 
pval <- cor.test.p(df); 

私たちは、すべての行にボンフェローニ補正を適用することにより、複数の仮説検証のための正しい(または列、行列が対称であるため)、我々は完了です。 p.adjustは、p値のベクトルを引数として取ります。

# Multiple hypothesis-testing corrected p-values 
# Note: pval is a symmetric matrix, so it doesn't matter if we correct 
# by column or by row 
padj <- apply(pval, 2, p.adjust, method = "bonferroni"); 
padj; 
#V1 V2  V3 V4  V5  V6 V7 V8  V9 V10 
#V1 0 1 1.0000000 1 1.0000000 1.0000000 1 1 1.0000000 1 
#V2 1 0 1.0000000 1 1.0000000 1.0000000 1 1 1.0000000 1 
#V3 1 1 0.0000000 1 0.9569498 1.0000000 1 1 1.0000000 1 
#V4 1 1 1.0000000 0 1.0000000 1.0000000 1 1 1.0000000 1 
#V5 1 1 0.9569498 1 0.0000000 1.0000000 1 1 1.0000000 1 
#V6 1 1 1.0000000 1 1.0000000 0.0000000 1 1 0.5461443 1 
#V7 1 1 1.0000000 1 1.0000000 1.0000000 0 1 1.0000000 1 
#V8 1 1 1.0000000 1 1.0000000 1.0000000 1 0 1.0000000 1 
#V9 1 1 1.0000000 1 1.0000000 0.5461443 1 1 0.0000000 1 
#V10 1 1 1.0000000 1 1.0000000 1.0000000 1 1 1.0000000 0 
+0

私は 'qdapTools'パッケージをダウンロードした後にやります。そのコードを実際に実行していると動作しませんでした。私が 'corr.test(dataHW8_1log、method = "pearson"、adjust = "bonferroni")を実行したときに、「何か重要なことがある場合にastreikが表示されず、不正確な結果が得られます。私が 'cor.test(dataHW8_1 $ Tree_Height、dataHW8_1 $ DBA、method = "pearson")'をチェックしたところ、 'p-value = 5.853e-06'の結果が得られました。しかし、私が 'corr.test'が私に与えるテーブルを見れば、確率値は '0.00'です。私はまだ平均、sdのテーブルを得る方法を失っています。 –

+0

@ N.J 'cor.test'はp値を与えます!アスタリスク(恐ろしい)アスタリスク表記は、必要に応じて自分で申請する必要があるものです。統計的有意性は、p値(あなたのケースでは、r = 0のピアソンの積 - 瞬間相関係数に対応するヌルの条件付き)に関して表現されます。私の例は間違いなく機能し、一貫しています。私は 'qdapTools'を知らないので、それを手伝うことはできません。例えばソースコードを使用するだけで、ソースデータの平均とsdを得ることができます。 'sapply(your_dataframe、mean)'と 'sapply(your_dataframe、sd)'です。 –

関連する問題