2016-08-18 1 views
1

私はある種の質問に対する参加者のバイナリの回答のデータセットを分析しています。私はDecの結果にVar * Base_conがどのように影響するかをテストするためにglm関数を使用しています。フィット後、 "Var"因子が "Base_con"因子レベルのそれぞれの結果にどのように影響するかを比較しようとしています。 this vignetteに続いて、私は私が再現できると信じて、次の(失敗)アプローチ、(それが動作しない場合は私に知らせてください)を試してみました:インタラクションの用語を含むGLMの後の複数の比較

# load example dataset with relevant columns 
require(RCurl) 
my_csv = getURL("https://docs.google.com/spreadsheets/d/1sBVW7QbnfumeRmY1uEDdiDNJE7QfmCXH0wMmV2lZSH4/pub?gid=0&single=true&output=csv") 
eg_data = read.csv(textConnection(my_csv)) 
# set columns as factors because they are numerically coded 
eg_data$Base_con = as.factor(eg_data$Base_con) 
eg_data$Var = as.factor(eg_data$Var) 
eg_data$Dec = as.factor(eg_data$Dec) 

# GLM fit 
m1 = glm(Dec ~ Var * Base_con, data = eg_data, family = "binomial") 

# strategy for Tukey multiple comparisons 
require(multcomp) 
tmp = expand.grid(Base_con = unique(eg_data$Base_con), Var = unique(eg_data$Var)) 
X = model.matrix(~Base_con : Var, data = tmp) 
mc = glht(m1, linfct = X) 

の最後のコマンドの出力は次のとおりです。

Error in glht.matrix(m1, linfct = X) : 
    ‘ncol(linfct)’ is not equal to ‘length(coef(model))’ 

効果は、エラーメッセージによって報告された2つの素子の列の数と長さが異なっている:私は前進することができたすべてのある

> ncol(X) 
[1] 7 
> length(coef(m1)) 
[1] 6 

これまでのところ。何か案は?ありがとうございます。あなたの目的はBase_conVarの相互作用の効果を比較する場合は、より直接的な方法は、これらの相互作用のための新しい用語を作成するかもしれ

+0

注意。 'glm()'は基本パッケージにあります。参加者ごとに複数の回答がある場合は、おそらく 'glmer'を使う必要があります。 – gung

+0

@gung、そうです。それは基本的なglm関数です。今回は繰り返し対策がないので、 'glmer'を使う必要はありません。 @ Silverfish、再現可能な例を使用して編集しようとしますが、エラーを再現できるかどうかはわかりません。 –

+0

*あなた*が@Lea_Casiraghiというエラーを再現できない場合、どうすれば*できますか?エラーを再現しないで解決策を見つけることができますか? – gung

答えて

1

(ここで私はdにデータフレーム名を変更した):

d$BV <- interaction(d$Base_con, d$Var) 

は次にモデルに適合して比較を行います。 `lme4`パッケージは* *` glm`機能を持っていないことを

# GLM fit 
m1 <- glm(Dec ~ -1 + BV, data = d, family = "binomial") 

library(multcomp) 
summary(glht(m1, linfct = mcp(BV = "Tukey"))) 

    # Simultaneous Tests for General Linear Hypotheses 

    # Multiple Comparisons of Means: Tukey Contrasts 


    # Fit: glm(formula = Dec ~ -1 + BV, family = "binomial", data = d) 

    # Linear Hypotheses: 
    #    Estimate Std. Error z value Pr(>|z|)  
    # 2.0 - 1.0 == 0 -1.7988  0.4632 -3.883 0.00133 ** 
    # 3.0 - 1.0 == 0 -4.9702  0.4846 -10.255 < 0.001 *** 
    # 1.1 - 1.0 == 0 -1.6596  0.4523 -3.669 0.00308 ** 
    # 2.1 - 1.0 == 0 -3.0593  0.4392 -6.965 < 0.001 *** 
    # 3.1 - 1.0 == 0 -5.3893  0.4759 -11.325 < 0.001 *** 
    # 3.0 - 2.0 == 0 -3.1714  0.3190 -9.941 < 0.001 *** 
    # 1.1 - 2.0 == 0 0.1392  0.2673 0.521 0.99498  
    # 2.1 - 2.0 == 0 -1.2605  0.2446 -5.154 < 0.001 *** 
    # 3.1 - 2.0 == 0 -3.5905  0.3055 -11.751 < 0.001 *** 
    # 1.1 - 3.0 == 0 3.3106  0.3029 10.930 < 0.001 *** 
    # 2.1 - 3.0 == 0 1.9109  0.2830 6.751 < 0.001 *** 
    # 3.1 - 3.0 == 0 -0.4191  0.3371 -1.243 0.80488  
    # 2.1 - 1.1 == 0 -1.3997  0.2231 -6.273 < 0.001 *** 
    # 3.1 - 1.1 == 0 -3.7297  0.2887 -12.920 < 0.001 *** 
    # 3.1 - 2.1 == 0 -2.3300  0.2678 -8.702 < 0.001 *** 
    # --- 
    # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    # (Adjusted p values reported -- single-step method) 
+0

私はその結果を再現することができましたが、私はそれが何を意味するのか、それが何を比較しているのか分かりません。 –

+0

6つのVarとBase_conの組み合わせが必要であり、合計15の可能な比較と、Base_conレベル内のVar値の3つの比較が必要です。私はここで何がテストされているのか、それとも私の例ではうまくいかないのかは分かりません。 –

+0

ここでの比較は、例えば、 'Var0:Base_con1 == 0'、つまり0との差のテスト。エフェクトを比較したい場合は、' X'を別に設定する必要があります。数分で何かを置くことができるかどうかがわかります。 –

関連する問題