ElMaestro ★★★ Denmark, 20090520 20:35 (edited by ElMaestro on 20090521 10:29) Posting: # 3719 Views: 9,368 

Dear all, I have two questions about R and BE and the model specification of Chow & Liu. In principle, I could ask these questions in the more dedicated R fora, however, since Bebac users are friendlier and sexier and probably are able to provide answers that relate directly to BE, I prefer to ask them here. For illustration purposes I have invented some arbitrary data, in R syntax:
lnY = c(10,11,12,13,14,15,10,11,12,13,14,15,16,17) ## seems to be standard syntax in this forum...which corresponds to a 2,2,2BE study with 7 subjects 1. Chow & Liu's model specification is lnY = Mean + S(i,k) + P(j) + F(j,k) + C(j1,k) + E(i,j,k) The lm fit above is not exactly corresponding to this model, as lm will fit an intercept unless we request it specifically not to. The fitted intercept is not the mean. This makes no difference in terms of the anova/error term BUT it will have some consequence for the regressors. Luckily, when we calculate the 90% CI, we do not need the individual model constants for the two treatments, we actually only need the difference between the model coefficients for the two treatments [rule: ln(A)ln(B)=ln(A/B)], and that difference is preserved regardless of the magnitude of the intercept. We can get the difference from summary(Snot) , where we read out how different Trt2 is from Trt1. But let's assume we want to do the fit as Chow & Liu want us to. Which lm syntax works? I cannot do it with lm or glm . The only way I can do achieve the same is via the absolutely horrendous syntax:
Snot2 = lm(lnYmean(lnY)~Trt+Subj+Seq+Per1) or:lnY2 = lnYmean(lnY) and then Snot2 = lm(lnY2~Trt+Subj+Seq+Per1) The advantage of Snot2 is that we ask R to do a fit without the intercept (we have subtracted the intercept/mean on the lefthand side anyway). Now we can read out both treatment regressors directly from summary(Snot2) , not just the difference between them. Is there a smarter, direct way of doing the fit the Chow&Liu way with mean as intercept? 2. In principle the overall model matrix for Snot above will have 14 rows (14 observations) and 14 columns (1 intercept+7 subjects+2 periods+2 treatments+2 sequences). Then there are redundancies so the 14×14 matrix will not be full rank, making it difficult to do matrix boogie (like inversion etc, which can be really handy when fitting) on it. Some statistical packages therefore prefer to reduce the model matrices to full rank. Now let's look at R's model matrix for Snot:
model.matrix(Snot) On my computer, I get a matrix with 10 columns. Obviously, the column called "Seq2" can be constructed from the three columns "Subj5", "Subj6" and "Subj7". Thus the muddle matrix is not full rank. It has rank=9 [ qr(model.matrix(Snot)) ]. So, why on earth will R reduce a model matrix into something that does not have full rank? What is the computational advantage of doing the muddle fit on this matrix? Many thanks for any insight you can provide EM.  Edit1: Formatted using BBCode. [Helmut] Edit2: Originally got 15 when I added 1,7,2,2 and 2. Corrected above. Morale: Never trust any calculation done by ElMaestro. Edit3: Welcome to the club! [Helmut] 
yjlee168 ★★ Kaohsiung, Taiwan, 20090521 13:35 @ ElMaestro Posting: # 3721 Views: 8,415 

dear Elmaestro, So far we still cannot figure out why the intercept is not the mean of lnY. However, you put Subj as one of regressors in your lm() functions. In Chow & Liu's model, the S(i,k) is the intersubject residuals of Subj i on Sequence k, not Subj itself. Could this be the reason that these two models might be different? What we use in bear with your example is like the following: Snot < lm(lnY~Trt+Subj:Seq+Per+Seq) for a 2x2x2 ABE. It works o.k..» Snot = lm(lnY~Trt+Subj+Seq+Per) ## seems to be standard syntax in » 1. Chow & Liu's model specification is » lnY = Mean + S(i,k) + P(j) + F(j,k) + C(j1,k) + E(i,j,k) — All the best, Yungjin Lee bear v2.8.4: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
ElMaestro ★★★ Denmark, 20090521 14:30 @ yjlee168 Posting: # 3722 Views: 8,380 

Hi Bears » Snot < lm(lnY~Trt+Subj:Seq+Per+Seq) for a 2x2x2 ABE. It» works o.k.. I guess you are right, but I think this is fishy. I could be wrong, but the colon as far as I understand from the manual is used for interactions?! Try and look at the model matrix with your syntax, where you will see columns filled with zeros. For example, with these data you will see a column called "Subj2:Seq2" which consists of all zeros, simply because Subj2 is assgined to Seq1 (and therefore the column called "Subj2:Seq1" has a few ones in it). A zero clumn is meaningless as it conveys no info to the fitting. In my naïvity I just thought Subj%in%Seq would be better, and I just tried it, but that seems to result in the same phenomenon anyway. In both cases if you compare e.g. the column "Subj2:Seq1" with the "Subj2" column of the original fit, they are identical (and so forth). Therefore, in this case I think the explicit declaration of a nested factor in the model is not necessary.EM. 
ElMaestro ★★★ Denmark, 20090522 14:39 @ ElMaestro Posting: # 3728 Views: 8,465 

Dear ElMaestro, I was lying awake all night trying to figure this one out. It appears that options(contrasts = c("contr.helmert", "contr.poly")) does the trick, then the model fits the mean as intercept. I got the inspiration for this here. Let it be said, I am not totally comfortable with contrasts (or with BE or with R or ...) but if the contrasting principle of summing to zero applies here then I would guess that the coefficient for Trt1 is Trt2; the summary readout is consistent with this.Aaaah, life is one long bloody learning process. HS, I guess Karl Popper or Kermit the Frog or some other fancy philosopher has a quote that fits on me in this situation?! Best regards ElMaestro  "Roses are red, violets are blue, I'm a schizophrenic and so am I" (attributed to Oscar Levant) 
yjlee168 ★★ Kaohsiung, Taiwan, 20090523 20:06 @ ElMaestro Posting: # 3741 Views: 8,367 

Dear ElMaestro, I don't quite understand the link that you mentioned in this post. So it's called "contrast coding"? And what you do mean about '...principle of summing to zero...'? We have already known that SPSS, SAS, and now R have different coding policy for dummy variables. Thanks, EM. » It appears that options(contrasts = c("contr.helmert", » "contr.poly")) does the trick, then the model fits the mean as» intercept. I got the inspiration for this here. Let it be said, I am » not totally comfortable with contrasts (or with BE or with R or ...) — All the best, Yungjin Lee bear v2.8.4: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
Aceto81 ★ Belgium, 20090526 10:28 @ ElMaestro Posting: # 3760 Views: 8,166 

There is another way: first substract the intercept, then add an offset, which is always fit with coefficient 1: If you take the mean as an offset, you get an error: "lengths differ", so if you use ave, you get: lm(formula = lnY ~ Trt + Subj + Seq + Per  1 + offset(ave(lnY))) Ace 
yjlee168 ★★ Kaohsiung, Taiwan, 20090523 19:30 @ ElMaestro Posting: # 3739 Views: 8,269 

Dear ElMaestro, We have tested R codes of your examples with SAS. And SAS comes out the correct answer with intercept equal to the mean. So we have been thinking what it is going wrong with R. Finally, we find that if we change the codes of dummy variables as follows: » Per = as.factor(c(1,1,1,1,1,1,1,2,2,2,2,2,2,2)) toPer = as.factor(c(1,1,1,1,1,1,1,0,0,0,0,0,0,0)) » Seq = as.factor(c(1,1,1,1,2,2,2,1,1,1,1,2,2,2)) toSeq = as.factor(c(1,1,1,1,0,0,0,1,1,1,1,0,0,0)) » Trt = as.factor(c(1,1,1,1,2,2,2,2,2,2,2,1,1,1)) toTrt = as.factor(c(1,1,1,1,0,0,0,0,0,0,0,1,1,1)) Then R and SAS have the same results. Apparently, R has somewhat different from SAS in dealing with lm() , or even with many others... There must be a term in statistics to describe this difference. We just don't know the term yet. However, it becomes more and more interesting now...— All the best, Yungjin Lee bear v2.8.4: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
ElMaestro ★★★ Denmark, 20090523 19:53 (edited by ElMaestro on 20090523 21:33) @ yjlee168 Posting: # 3740 Views: 8,300 

Dear Bears » We just don't know the term yet. "Contrast coding", is my guess. I think R uses a different method of contrasts in the model matrix. I don't have SAS so I have idea if this is correct. Try and look at SAS' model matrices yourself. Contrasts are the secret to how the fit is done; one model can be fit with different contasts, and the model coefficients depend on the contrast method. In the default setting R puts a one in the m.m. in a given column when the corresponding dataline contains this factor at the level indicated by the column. So, for example, in our example the columns corresponding to subjects will sum to two because we have two observations from each subject. Another way of contrasting is make them sum to zero. If we want to be sure we get the model is fit with the mean as intercept, I am now inclined (bear in mind I am just an amateur equipped with a brain having the size of a walnut) to think that we need the contrasts to sum to zero; all errors in any model by definition must sum to zero. In Chow&Liu's model, mean plus all coefficients then must sum to the mean, so all coefficient must sum to zero. EM. 
d_labes ★★★ Berlin, Germany, 20090525 11:47 @ ElMaestro Posting: # 3746 Views: 8,268 

Dear ElMaestro, dear bears, Let me throw my 2 cents into this discussion. What are we talking about? If I understand it, about reparameterization to solve the equations for the coefficients of the ANOVA model(s). Recall your very first lessons about ANOVA (as I had to do) to remember that the ANOVA model in case of categorical effects (factors) is overparameterized. To solve nevertheless constraints on the coefficients of the model have to imposed. Historical the first? and very common, covered in many statistical textbooks, is the "sum to zero" reparameterization, i.e. all coefficients for the various levels of a distinct factor sum to zero. This results in the intercept being the overall mean of the independent variable as ElMaestro noticed (using a variant of this, the so called "Helmert contrast", which I never had understood fully ). The correct specification in R for the simple "sum to zero" is options(contrasts = c("contr.sum", "contr.poly")) .Other common reparameterization is the "set to zero". This is done by simply setting the coefficient of one (reference) level of a factor to zero. R slang for that is options(contrasts = c("contr.treatment", "contr.poly")) .This is also the setting of "The power to know" with the solely difference that R / S and other statistical software sets the first level to zero and SAS the last. As always SAS is outstanding . This is the reason why Yungjin has noticed the same results in R and SAS after recoding the last level being the first. BTW you can obtain the same results without recoding if you use library(nlme) # needed for the contrast contr.SAS The "set to zero" reparameterization yields an intercept, that is not! in general the overall mean of the independent variable, as Yungjin's statement suggests: » We have tested R codes of your examples with SAS. And SAS comes out the » correct answer with intercept equal to the mean. [...] Therefore, Yungjin, check your SAS results. There are other reparameterizations out there. But do not ask me why and what they really mean and what they are good for. BTW: Don't ask me further, what "contr.poly" stands for in the R code above. I'm a very beginneR. BTW2: Dear ElMaestro. Intelligence is not a linear function of brain size nor brain weight as neuro science has long ago recognized. The German saying "Leave the thinking to the horses, they have a greater head" is definitely not true! BTW3: IMHO there is no definite advantage or disadvantage of the reparameterization with respect to the BE test, as you can see if you calculate the difference and CI of factor levels. Only the interpretation of the "raw" coefficients of the ANOVA model are changed. — Regards, Detlew 
yjlee168 ★★ Kaohsiung, Taiwan, 20090525 14:17 @ d_labes Posting: # 3747 Views: 8,297 

Dear D. Labes and ElMaestro, Thank you for your reminding. I just checked mt SAS output as the follows.
The SAS System And output from R: Call: Yes, you're correct about this. The results are the same from R and SAS. However, it's not the overall mean for lnY. Sorry about this. » » We have tested R codes of your examples with SAS. And SAS comes out the » » correct answer with intercept equal to the mean. [...]» » » Therefore, Yungjin, check your SAS results. — All the best, Yungjin Lee bear v2.8.4: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
ElMaestro ★★★ Denmark, 20090525 20:29 @ d_labes Posting: # 3752 Views: 8,274 

Dear dlabes, » BTW2: Dear ElMaestro. Intelligence is not a linear function of brain size nor brain weight as neuro science has long ago recognized. Oh. That means I cannot even claim intellectual noninferiority to e.g. a skunk . Thank you. EM. 
d_labes ★★★ Berlin, Germany, 20090526 08:01 @ ElMaestro Posting: # 3756 Views: 8,241 

Dear ElMaestro, » Oh. That means I cannot even claim intellectual noninferiority to e.g. a skunk . Thank you. Think positive! It means a brain with size of a walnut can have intellectual superiority over an potwhale . — Regards, Detlew 