अवलोकन द्वारा मैं हल्के निश्चित प्रभाव कैसे निकालें?

Question

मैं एक lme वस्तु, कुछ दोहराया उपायों से निर्मित है पोषक तत्व सेवन डेटा (RespondentID प्रति दो 24 घंटे का सेवन अवधि):

Male.lme2 <- lmer(BoxCoxXY ~ -1 + AgeFactor + IntakeDay + (1|RespondentID), 
    data = Male.Data, 
    weights = SampleWeight) 

और मैं सफलतापूर्वक ranef(Male.lme1) का उपयोग कर RespondentID द्वारा यादृच्छिक प्रभाव प्राप्त कर सकते हैं। मैं RespondentID द्वारा निश्चित प्रभावों का परिणाम भी एकत्र करना चाहता हूं। जैसा कि मैंने नीचे दिखाया है, coef(Male.lme1) बिल्कुल वही प्रदान नहीं करता है जो मुझे चाहिए।

> summary(Male.lme1) 
Linear mixed model fit by REML 
Formula: BoxCoxXY ~ AgeFactor + IntakeDay + (1 | RespondentID) 
    Data: Male.Data 
    AIC BIC logLik deviance REMLdev 
    9994 10039 -4990  9952 9980 
Random effects: 
Groups  Name  Variance Std.Dev. 
RespondentID (Intercept) 0.19408 0.44055 
Residual     0.37491 0.61230 
Number of obs: 4498, groups: RespondentID, 2249 

Fixed effects: 
        Estimate Std. Error t value 
(Intercept)   13.98016 0.03405 410.6 
AgeFactor4to8  0.50572 0.04084 12.4 
AgeFactor9to13  0.94329 0.04159 22.7 
AgeFactor14to18  1.30654 0.04312 30.3 
IntakeDayDay2Intake -0.13871 0.01809 -7.7 

Correlation of Fixed Effects: 
      (Intr) AgFc48 AgF913 AF1418 
AgeFactr4t8 -0.775      
AgeFctr9t13 -0.761 0.634    
AgFctr14t18 -0.734 0.612 0.601  
IntkDyDy2In -0.266 0.000 0.000 0.000 

मैं अपने डेटा को फिट परिणाम संलग्न है, head(Male.Data) से पता चलता

NutrientID RespondentID Gender Age SampleWeight IntakeDay IntakeAmt AgeFactor BoxCoxXY lmefits 
2   267  100020  1 12 0.4952835 Day1Intake 12145.852  9to13 15.61196 15.22633 
7   267  100419  1 14 0.3632839 Day1Intake 9591.953 14to18 15.01444 15.31373 
8   267  100459  1 11 0.4952835 Day1Intake 7838.713  9to13 14.51458 15.00062 
12  267  101138  1 15 1.3258785 Day1Intake 11113.266 14to18 15.38541 15.75337 
14  267  101214  1 6 2.1198688 Day1Intake 7150.133  4to8 14.29022 14.32658 
18  267  101389  1 5 2.1198688 Day1Intake 5091.528  4to8 13.47928 14.58117 

coef(Male.lme1) से लाइनों के पहले दो हैं:

$RespondentID 
     (Intercept) AgeFactor4to8 AgeFactor9to13 AgeFactor14to18 IntakeDayDay2Intake 
100020 14.28304  0.5057221  0.9432941  1.306542   -0.1387098 
100419 14.00719  0.5057221  0.9432941  1.306542   -0.1387098 
100459 14.05732  0.5057221  0.9432941  1.306542   -0.1387098 
101138 14.44682  0.5057221  0.9432941  1.306542   -0.1387098 
101214 13.82086  0.5057221  0.9432941  1.306542   -0.1387098 
101389 14.07545  0.5057221  0.9432941  1.306542   -0.1387098 

प्रदर्शित करने के लिए कैसे coef परिणाम से संबंधित Male.Data में फिट अनुमान (जिन्हें Male.Data$lmefits <- fitted(Male.lme1) का उपयोग करके पकड़ा गया था, पहले RespondentID के लिए, जिसकी आयु है फैक्टर स्तर 9-13: - सज्जित मूल्य 15.22633 है, जो बराबर होती है - coeffs से - (Intercept) + (AgeFactor9-13) = 14.28304 + 0.9432941

वहाँ एक चतुर आदेश मुझे इस्तेमाल के लिए है कि क्या करेंगे मैं अपने आप चाहते हैं, जो निश्चित प्रभाव को निकालने के लिए है प्रत्येक विषय के लिए अनुमान, या क्या मुझे if की श्रृंखला का सामना करना पड़ रहा है, जिसमें इंटरसेप्ट से यादृच्छिक प्रभाव योगदान को कम करने के बाद, सही निश्चित प्रभाव अनुमान प्राप्त करने के लिए प्रत्येक विषय में सही आयुफैक्टर स्तर लागू करने का प्रयास किया गया है?

अद्यतन, क्षमा, मैं आउटपुट पर कटौती करने की कोशिश कर रहा था और str() के बारे में भूल गया था। आउटपुट है:

>str(Male.Data) 
'data.frame': 4498 obs. of 11 variables: 
$ NutrientID : int 267 267 267 267 267 267 267 267 267 267 ... 
$ RespondentID: Factor w/ 2249 levels "100020","100419",..: 1 2 3 4 5 6 7 8 9 10 ... 
$ Gender  : int 1 1 1 1 1 1 1 1 1 1 ... 
$ Age   : int 12 14 11 15 6 5 10 2 2 9 ... 
$ BodyWeight : num 51.6 46.3 46.1 63.2 28.4 18 38.2 14.4 14.6 32.1 ... 
$ SampleWeight: num 0.495 0.363 0.495 1.326 2.12 ... 
$ IntakeDay : Factor w/ 2 levels "Day1Intake","Day2Intake": 1 1 1 1 1 1 1 1 1 1 ... 
$ IntakeAmt : num 12146 9592 7839 11113 7150 ... 
$ AgeFactor : Factor w/ 4 levels "1to3","4to8",..: 3 4 3 4 2 2 3 1 1 3 ... 
$ BoxCoxXY : num 15.6 15 14.5 15.4 14.3 ... 
$ lmefits  : num 15.2 15.3 15 15.8 14.3 ... 

वजन और लिंग इस्तेमाल किया जा रहा नहीं कर रहे हैं (यह पुरुषों डेटा है, इसलिए सभी लिंग मान ही कर रहे हैं) और NutrientID इसी तरह के डेटा के लिए तय हो गई है।

मैं भयानक ifelse कथन कर रहा हूं जो मैंने पोस्ट किया है, इसलिए तुरंत आपके सुझाव का प्रयास करेंगे। :)

अद्यतन 2: यह मेरे वर्तमान डेटा के साथ पूरी तरह से काम करता है और इसके लिए टिप्पणी में अतिरिक्त सहायता के लिए डीडब्ल्यूएन के लिए धन्यवाद, नए डेटा के लिए भविष्य के सबूत होना चाहिए। :)

AgeLevels <- length(unique(Male.Data$AgeFactor)) 
Temp <- as.data.frame(fixef(Male.lme1)['(Intercept)'] + 
c(0,fixef(Male.lme1)[2:AgeLevels])[ 
     match(Male.Data$AgeFactor, c("1to3", "4to8", "9to13","14to18", "19to30","31to50","51to70","71Plus"))] + 
c(0,fixef(Male.lme1)[(AgeLevels+1)])[ 
     match(Male.Data$IntakeDay, c("Day1Intake","Day2Intake"))]) 
names(Temp) <- c("FxdEffct") 

Answer 1

यह यद्यपि आप वास्तव में हमें दिया है चाहिए str के परिणाम (Male.Data) कुछ इस तरह (होने जा रहा है क्योंकि मॉडल उत्पादन हमें आधारभूत मूल्यों के लिए कारक स्तरों नहीं बताता है :)

#First look at the coefficients 
fixef(Male.lme2) 

#Then do the calculations 
fixef(Male.lme2)[`(Intercept)`] + 
c(0,fixef(Male.lme2)[2:4])[ 
      match(Male.Data$AgeFactor, c("1to3", "4to8", "9to13","14to18"))] + 
c(0,fixef(Male.lme2)[5])[ 
      match(Male.Data$IntakeDay, c("Day1Intake","Day2Intake"))] 

आप मूल रूप से एक match समारोह के माध्यम से मूल डेटा चल रहे हैं सही गुणांक (रों) अवरोधन में जोड़ने के लिए लेने के लिए ... 0 हो जाएगा जो करता है, तो डेटा कारक के आधार स्तर (जिसका है वर्तनी मैं अनुमान लगा रहा हूं।)

संपादित करें: मैंने अभी देखा है कि आपने फॉर्मूला में "-1" लगाया है, इसलिए शायद आपके सभी आयुफैक्टर शब्द आउटपुट में सूचीबद्ध हैं और आप गुणांक वेक्टर में 0 और मिलान तालिका में आविष्कार किए गए आयुफैक्टर स्तर को बता सकते हैं वेक्टर।

Answer 2

नीचे यह है कि lme4 -package में व्यक्तियों के निश्चित प्रभाव और यादृच्छिक प्रभाव घटकों को निकालने के लिए मुझे हमेशा सबसे आसान लगता है। यह वास्तव में प्रत्येक अवलोकन के लिए उपयुक्त फिट निष्कर्ष निकालता है। मान लिया जाये कि हम फार्म की एक मिश्रित प्रभाव मॉडल है:

y = Xb + Zu + e 

जहां Xb तय इफेक्ट होते हैं और जेड यू यादृच्छिक प्रभाव होते हैं, हम घटकों को निकालने कर सकते हैं (एक उदाहरण के रूप lme4 के sleepstudy उपयोग करते हुए):

library(lme4) 
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) 

# Xb 
fix <- getME(fm1,'X') %*% fixef(fm1) 
# Zu 
ran <- t(as.matrix(getME(fm1,'Zt'))) %*% unlist(ranef(fm1)) 
# Xb + Zu 
fixran <- fix + ran 

मुझे पता है कि यह रैखिक मिश्रित प्रभाव मॉडल से घटकों को निकालने के लिए एक सामान्यीकृत दृष्टिकोण के रूप में कार्य करता है। गैर-रैखिक मॉडल के लिए, मॉडल मैट्रिक्स एक्स में दोहराव होता है और आपको उपरोक्त कोड को थोड़ा सा बनाना पड़ सकता है। यहाँ कुछ मान्यता उत्पादन के साथ ही एक जाली का उपयोग कर दृश्य है:

> head(cbind(fix, ran, fixran, fitted(fm1))) 
     [,1]  [,2]  [,3]  [,4] 
[1,] 251.4051 2.257187 253.6623 253.6623 
[2,] 261.8724 11.456439 273.3288 273.3288 
[3,] 272.3397 20.655691 292.9954 292.9954 
[4,] 282.8070 29.854944 312.6619 312.6619 
[5,] 293.2742 39.054196 332.3284 332.3284 
[6,] 303.7415 48.253449 351.9950 351.9950 

# Xb + Zu 
> all(round((fixran),6) == round(fitted(fm1),6)) 
[1] TRUE 

# e = y - (Xb + Zu) 
> all(round(resid(fm1),6) == round(sleepstudy[,"Reaction"]-(fixran),6)) 
[1] TRUE 

nobs <- 10 # 10 observations per subject 
legend = list(text=list(c("y", "Xb + Zu", "Xb")), lines = list(col=c("blue", "red", "black"), pch=c(1,1,1), lwd=c(1,1,1), type=c("b","b","b"))) 
require(lattice) 
xyplot(
    Reaction ~ Days | Subject, data = sleepstudy, 
    panel = function(x, y, ...){ 
     panel.points(x, y, type='b', col='blue') 
     panel.points(x, fix[(1+nobs*(panel.number()-1)):(nobs*(panel.number()))], type='b', col='black') 
     panel.points(x, fixran[(1+nobs*(panel.number()-1)):(nobs*(panel.number()))], type='b', col='red') 
    }, 
    key = legend 
)