lmer의 절편 선택

Question

lmerTest의 lmer 모델을 사용하여 크기가 유전자 발현과 유의미한 상관 관계가 있는지 알아보고 어떤 경우에는 어떤 특정 유전자가 크기와 관련이 있는지 알아 봅니다 ('여성'및 '새장'도 고려함). '랜덤 효과) : 그것은, 알파벳에서 가장 높은이기 때문에 요약 출력에서 ​​

lmer(Expression ~ size*genes + (1|female) + (1|cage), data = df) 

, 내 유전자 중 하나가 (절편으로까지 사용되는' 'ctsk). 주변을 읽은 후에는 다른 모든 것을 비교하기 위해 가장 높은 (또는 가장 낮은) 발현 유전자를 내 도청 물로 선택하는 것이 좋습니다. 이 경우 유전자 '스타'가 가장 높게 표현되었습니다. 내 데이터의 레벨을 다시 조정하고 'star'를 절편으로 사용하여 모델을 다시 실행 한 후에는 anova() 출력이 동일하지만 summary() 출력에서 ​​다른 모든 경사가 중요합니다.

내 질문은 :

1. 은 절편으로 사용 내 유전자 중 하나가되지 수 있습니까? 그것이 가능하지 않다면, 어떤 유전자를 절편으로 선택해야하는지 어떻게 알 수 있습니까?
2. 슬로프가 0과 다른지 테스트 할 수 있습니까? 아마도 이것이 내 모델 (예 : '0 + 크기 * 유전자')에서 절편을 지정하지 않을 것입니다.
3. 절편을 모든 경사면의 평균으로 가질 수 있습니까?

다음으로 lsmeans를 사용하여 기울기가 서로 다른지 여부를 결정합니다. 그렇지 않으면-희망이 괜찮

df <- structure(list(size = c(13.458, 13.916, 13.356, 13.84, 14.15, 
          16.4, 15.528, 13.916, 13.458, 13.285, 15.415, 14.181, 13.367, 
          13.356, 13.947, 14.615, 15.804, 15.528, 16.811, 14.677, 13.2, 
          17.57, 13.947, 14.15, 16.833, 13.2, 17.254, 16.4, 14.181, 13.367, 
          14.294, 13.84, 16.833, 17.083, 15.847, 13.399, 14.15, 15.47, 
          13.356, 14.615, 15.415, 15.596, 15.847, 16.833, 13.285, 15.47, 
          15.596, 14.181, 13.356, 14.294, 15.415, 15.363, 15.4, 12.851, 
          17.254, 13.285, 17.57, 14.7, 17.57, 13.947, 16.811, 15.4, 13.399, 
          14.22, 13.285, 14.344, 17.083, 15.363, 14.677, 15.945), female = structure(c(7L, 
                             12L, 7L, 11L, 12L, 9L, 6L, 12L, 7L, 7L, 6L, 12L, 8L, 7L, 7L, 
                             11L, 9L, 6L, 10L, 11L, 8L, 10L, 7L, 12L, 10L, 8L, 10L, 9L, 12L, 
                             8L, 12L, 11L, 10L, 10L, 9L, 8L, 12L, 6L, 7L, 11L, 6L, 9L, 9L, 
                             10L, 7L, 6L, 9L, 12L, 7L, 12L, 6L, 6L, 6L, 8L, 10L, 7L, 10L, 
                             11L, 10L, 7L, 10L, 6L, 8L, 11L, 7L, 6L, 10L, 6L, 11L, 9L), .Label = c("2", 
                                              "3", "6", "10", "11", "16", "18", "24", "25", "28", "30", "31", 
                                              "116", "119", "128", "135", "150", "180", "182", "184", "191", 
                                              "194", "308", "311", "313", "315", "320", "321", "322", "324", 
                                              "325", "329", "339", "342"), class = "factor"), Expression = c(1.10620339407889, 
                                                              1.06152707257767, 2.03000185674761, 1.92971750056866, 1.30833983462599, 
                                                              1.02760836165184, 0.960969703469363, 1.54706275342441, 0.314774666283256, 
                                                              2.63330873720495, 0.895123048920455, 0.917716470037954, 1.3178821021651, 
                                                              1.57879156856332, 0.633429011784367, 1.12641940390116, 1.0117475796626, 
                                                              0.687813581350802, 0.923485880847423, 2.98926377892241, 0.547685277701021, 
                                                              0.967691178046748, 2.04562285257417, 1.09072264997544, 1.57682235413366, 
                                                              0.967061529758701, 0.941995966023426, 0.299517719292817, 1.8654758451133, 
                                                              0.651369936708288, 1, 1.04407979584122, 0.799275069735012, 1.007255409328, 
                                                              0.428129727802404, 0.93927930755046, 0.987394257033815, 0.965050972503591, 
                                                              2.06719308587322, 1.63846508102874, 0.997380526962644, 0.60270197593643, 
                                                              2.78682867333149, 0.552922632281237, 3.06702198884562, 0.890708510580522, 
                                                              1.15168812515828, 0.929205084743164, 2.27254101826041, 1, 0.958147442333527, 
                                                              1.05924173014089, 0.984356852670054, 0.623630720815415, 0.796864961771971, 
                                                              2.4679841984147, 1.07248904053777, 1.79630829771291, 0.929642913565982, 
                                                              0.296954006040077, 2.25741254504115, 1.17188536743493, 0.849778293699644, 
                                                              2.32679163466857, 0.598119006609413, 0.975660099975423, 1.01494421228949, 
                                                              1.14007557533352, 2.03638316428189, 0.777347547080068), cage = structure(c(64L, 
                                                                                 49L, 56L, 66L, 68L, 48L, 53L, 49L, 64L, 56L, 55L, 68L, 80L, 56L, 
                                                                                 64L, 75L, 69L, 53L, 59L, 66L, 63L, 59L, 64L, 68L, 59L, 63L, 50L, 
                                                                                 48L, 68L, 80L, 49L, 66L, 59L, 50L, 48L, 63L, 68L, 62L, 56L, 75L, 
                                                                                 55L, 81L, 48L, 59L, 56L, 62L, 81L, 68L, 56L, 49L, 55L, 62L, 55L, 
                                                                                 63L, 50L, 56L, 59L, 75L, 59L, 64L, 59L, 55L, 63L, 66L, 56L, 53L, 
                                                                                 50L, 62L, 66L, 81L), .Label = c("023", "024", "041", "042", "043", 
                                                                                         "044", "044 bis", "045", "046", "047", "049", "051", "053", "058", 
                                                                                         "060", "061", "068", "070", "071", "111", "112", "113", "123", 
                                                                                         "126", "128", "14", "15", "23 bis", "24", "39", "41", "42", "44", 
                                                                                         "46 bis", "47", "49", "51", "53", "58", "60", "61", "67", "68", 
                                                                                         "70", "75", "76", "9", "D520", "D521", "D522", "D526", "D526bis", 
                                                                                         "D533", "D535", "D539", "D544", "D545", "D545bis", "D546", "D561", 
                                                                                         "D561bis", "D564", "D570", "D581", "D584", "D586", "L611", "L616", 
                                                                                         "L633", "L634", "L635", "L635bis", "L637", "L659", "L673", "L676", 
                                                                                         "L686", "L717", "L718", "L720", "L725", "L727", "L727bis"), class = "factor"), 
       genes = c("igf1", "gr", "ctsk", "ets2", "ctsk", "mtor", "igf1", 
          "sgk1", "sgk1", "ghr1", "ghr1", "gr", "ctsk", "ets2", "timp2", 
          "timp2", "ets2", "rictor", "sparc", "mmp9", "gr", "sparc", 
          "mmp2", "ghr1", "mmp9", "sparc", "mmp2", "timp2", "star", 
          "sgk1", "mmp2", "gr", "mmp2", "rictor", "timp2", "mmp2", 
          "mmp2", "mmp2", "mmp2", "rictor", "mtor", "ghr1", "star", 
          "igf1", "mmp9", "igf1", "igf2", "rictor", "rictor", "mmp9", 
          "ets2", "ctsk", "mtor", "ghr1", "mtor", "ets2", "ets2", "igf2", 
          "igf1", "sgk1", "sgk1", "ghr1", "sgk1", "igf2", "star", "mtor", 
          "igf2", "ghr1", "mmp2", "rictor")), .Names = c("size", "female", 
                      "Expression", "cage", "genes"), row.names = c(1684L, 2674L, 10350L, 
                                 11338L, 10379L, 4586L, 1679L, 3637L, 3610L, 5537L, 5530L, 2676L, 
                                 10355L, 11313L, 8422L, 8450L, 11322L, 6494L, 9406L, 13262L, 2653L, 
                                 9407L, 12274L, 5564L, 13256L, 9394L, 12294L, 8438L, 750L, 3614L, 
                                 12303L, 2671L, 12293L, 6513L, 8437L, 12284L, 12305L, 12267L, 
                                 12276L, 6524L, 4567L, 5545L, 733L, 1700L, 13241L, 1674L, 7471L, 
                                 6528L, 6498L, 13266L, 11308L, 10347L, 4566L, 5541L, 4590L, 11315L, 
                                 11333L, 7482L, 1703L, 3607L, 3628L, 5529L, 3617L, 7483L, 722L, 
                                 4565L, 7476L, 5532L, 12299L, 6510L), class = "data.frame") 

genes <- as.factor(df$genes) 
library(lmerTest) 

fit1 <- lmer(Expression ~ size * genes +(1|female) + (1|cage), data = df) 
anova(fit1) 
summary(fit1) # uses the gene 'ctsk' for intercept, so re-level to see what happens if I re-order based on highest value (sgk1): 
df$genes <- relevel(genes, "star") 

# re-fit the model with 'star' as the intercept: 
fit1 <- lmer(Expression ~ size * genes +(1|female) + (1|cage), data = df) 
anova(fit1) # no difference here 
summary(fit1) # lots of difference 

내 샘플 데이터 꽤 긴 모델이 실행되지 않을 것이기 때문입니다 :

는 여기에 몇 가지 재현 코드입니다!

Answer 1

적합 모델의 계수를 해석하는 것이 가능하지만 가장 유익하고 생산적인 방법은 아닙니다. 대신 기본적으로 사용되는 명암 대비 방법을 사용하여 모델을 맞추고 적절한 사후 분석으로 후속 조치를 취하십시오. 이를 위해

, 내가 미래의 모든 발전이 일어날 것이다 lsmeans의 연속 인 emmeans (추정 한계 수단) 패키지를 사용하는 것이 좋습니다. 패키지에는 몇 개의 비 네트가 있으며 상황에 가장 적합한 것은 vignette("interactions")입니다. here - 특히 공변량과의 상호 작용 섹션을 볼 수 있습니다.

간단히 말하자면, 가로 채기를 비교하는 것은 보간법 인 size = 0의 예측이기 때문에 매우 오도 할 수 있습니다. 또한 질문에서 제안하는 것처럼 여기에서 요점은 아마도 요격을 요격보다 더 비교하는 것입니다. 이를 위해 emtrends() 함수가 있습니다 (또는 원하는 경우 별칭 lstrends()).

어떤 일이 벌어지고 있는지 시각화 할 수 있도록 모델 예측 그래프를 표시하는 것이 좋습니다.

library(emmeans) 
emmip(fit1, gene ~ size, at = list(size = range(df$size)))