I am in interested in how do effect coding in R. I know that someone else has asked this question (i.e How to do regression with effect coding instead of dummy coding in R?). Here is the lm() model on its own:
Fixed.Full<-lm(sqrtClay_Tot ~ Physio.Code+ROCKTYPE1, data=PedonsTx_Reread,
contrasts = list(Factor.0.Physio = contr.sum))
> summary(Fixed.Full)
Call:
lm(formula = sqrtClay_Tot ~ Physio.Code + ROCKTYPE1, data = PedonsTx_Reread)
Residuals:
Min 1Q Median 3Q Max
-4.5719 -1.1543 0.0102 1.1366 6.0146
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.18961 0.54194 9.576 < 2e-16 ***
Physio.Code2 -0.90588 0.15752 -5.751 1.11e-08 ***
Physio.Code3 0.70532 0.21141 3.336 0.000874 ***
Physio.Code4 0.56462 0.31295 1.804 0.071444 .
Physio.Code7 -0.04656 0.26410 -0.176 0.860075
Physio.Code8 -0.41242 0.24330 -1.695 0.090294 .
Physio.Code9 -0.81833 0.32048 -2.553 0.010782 *
Physio.Code11 -0.04619 0.16420 -0.281 0.778515
ROCKTYPE1chert 1.43390 1.80191 0.796 0.426319
ROCKTYPE1clay or mud 0.35031 0.53861 0.650 0.515562
ROCKTYPE1evaporite -0.95719 0.65049 -1.471 0.141408
ROCKTYPE1fine-grained mixed clastic -1.02193 0.59167 -1.727 0.084377 .
ROCKTYPE1gravel -0.16708 0.56139 -0.298 0.766041
ROCKTYPE1limestone 0.36667 0.55323 0.663 0.507585
ROCKTYPE1mudstone 0.18934 0.58553 0.323 0.746476
ROCKTYPE1rhyolite -0.25743 0.83684 -0.308 0.758426
ROCKTYPE1sand -0.92281 0.53193 -1.735 0.083014 .
ROCKTYPE1sandstone -1.57669 0.57599 -2.737 0.006281 **
ROCKTYPE1serpentinite -0.18647 1.84384 -0.101 0.919464
ROCKTYPE1shale 0.34408 0.62577 0.550 0.582518
ROCKTYPE1silt -0.82415 0.65554 -1.257 0.208912
ROCKTYPE1siltstone -1.21094 0.79293 -1.527 0.126970
ROCKTYPE1terrace -0.64082 0.60119 -1.066 0.286663
ROCKTYPE1tuff 1.10034 1.31108 0.839 0.401482
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.748 on 1261 degrees of freedom
Multiple R-squared: 0.2072, Adjusted R-squared: 0.1928
F-statistic: 14.33 on 23 and 1261 DF, p-value: < 2.2e-16
Now here is the model using effecting coding within lm():
> Fixed.Full<-lm(sqrtClay_Tot ~ Physio.Code+ROCKTYPE1, data=PedonsTx_Reread,
+ contrasts = list(Physio.Code = contr.sum))
> summary(Fixed.Full)
Call:
lm(formula = sqrtClay_Tot ~ Physio.Code + ROCKTYPE1, data = PedonsTx_Reread,
contrasts = list(Physio.Code = contr.sum))
Residuals:
Min 1Q Median 3Q Max
-4.5719 -1.1543 0.0102 1.1366 6.0146
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.06968 0.50933 9.954 < 2e-16 ***
Physio.Code1 0.11993 0.13291 0.902 0.36703
Physio.Code2 -0.78595 0.11689 -6.724 2.68e-11 ***
Physio.Code3 0.82525 0.15844 5.208 2.22e-07 ***
Physio.Code4 0.68455 0.24897 2.750 0.00605 **
Physio.Code5 0.07337 0.21084 0.348 0.72792
Physio.Code6 -0.29249 0.18905 -1.547 0.12208
Physio.Code7 -0.69840 0.26157 -2.670 0.00768 **
ROCKTYPE1chert 1.43390 1.80191 0.796 0.42632
ROCKTYPE1clay or mud 0.35031 0.53861 0.650 0.51556
ROCKTYPE1evaporite -0.95719 0.65049 -1.471 0.14141
ROCKTYPE1fine-grained mixed clastic -1.02193 0.59167 -1.727 0.08438 .
ROCKTYPE1gravel -0.16708 0.56139 -0.298 0.76604
ROCKTYPE1limestone 0.36667 0.55323 0.663 0.50759
ROCKTYPE1mudstone 0.18934 0.58553 0.323 0.74648
ROCKTYPE1rhyolite -0.25743 0.83684 -0.308 0.75843
ROCKTYPE1sand -0.92281 0.53193 -1.735 0.08301 .
ROCKTYPE1sandstone -1.57669 0.57599 -2.737 0.00628 **
ROCKTYPE1serpentinite -0.18647 1.84384 -0.101 0.91946
ROCKTYPE1shale 0.34408 0.62577 0.550 0.58252
ROCKTYPE1silt -0.82415 0.65554 -1.257 0.20891
ROCKTYPE1siltstone -1.21094 0.79293 -1.527 0.12697
ROCKTYPE1terrace -0.64082 0.60119 -1.066 0.28666
ROCKTYPE1tuff 1.10034 1.31108 0.839 0.40148
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.748 on 1261 degrees of freedom
Multiple R-squared: 0.2072, Adjusted R-squared: 0.1928
F-statistic: 14.33 on 23 and 1261 DF, p-value: < 2.2e-16
The code above is functional within the lm() function as seen above,and what I am trying to do in the code above is to predict clay content (square-root transformed) from the categorical variables 'Physio.Code' and 'ROCKTYPE1'. What I would like to know, however, is if there is any way to do effect coding using the gls() or glm() functions in R, where method="REML" or REML estimation can be used. I am interested in comparing the levels within 'Physio.Code' specifically to to the grand sample mean (i.e. sum contrasts).
The other question I had was how to interpret the numbers after the categorical variable and whether they correspond to the same levels. With respect to 'Physio.Code', in the two models above '1' was used as the baseline - the levels were originally '1,2,3,4,7,8,9,11" but when '1' was used as the baseline group, it switched to "2,3,4,7,8,9" - I wasn't sure if levels are the same (i.e. Physio.Code2 in model is the same as Physio.Code1 in the 2nd model). Thank you.
