0
$\begingroup$

I have my data that looks like this: data: enter image description here

or in R:

df <- read.table(header= TRUE, text=
                   "Location    treatment   block   chem1   chem2   chem3   yield
Monmouth    CC  1   1.600946    6212.631579 810.3024    5156
Monmouth    CC  2   1.474084    5641.403509 835.5242    5157
Monmouth    CC  3   1.433078    4392.280702 856.206 5551
Monmouth    CC  4   1.532492    7120.701754 759.8589    5156
Monmouth    CCS 1   1.198667    5305.964912 726.2298    4804
Monmouth    CCS 2   1.193193    4936.842105 792.6472    4720
Monmouth    CCS 3   1.214941    5383.157895 724.7165    4757
Monmouth    CCS 4   1.360981    6077.894737 699.3266    5021
Monmouth    CS  1   1.853056    5769.122807 500.9153    4826
Monmouth    CS  2   1.690117    5016.842105 634.8698    4807
Monmouth    CS  3   1.544357    5060.350877 637.9536    4475
Monmouth    CS  4   1.825132    5967.017544 648.8814    4596
Monmouth    SCS 1   1.695409    5576.842105 641.0357    4669
Monmouth    SCS 2   1.764123    5174.035088 623.3822    4588
Monmouth    SCS 3   1.903743    5655.438596 555.2834    4542
Monmouth    SCS 4   1.538684    5468.77193  520.8119    4592
Urbana  CC  1   0.845077    7933.333333 683.3528    5583
Urbana  CC  2   1.011463    7000    595.9173    5442
Urbana  CC  3   0.857032    6352.982456 635.4315    5693
Urbana  CC  4   0.989803    8153.684211 672.4234    5739
Urbana  CCS 1   0.859022    6077.894737 847.2944    5045
Urbana  CCS 2   0.919467    4939.649123 745.5665    4902
Urbana  CCS 3   1.01717 5002.807018 778.3548    5006
Urbana  CCS 4   0.861689    6489.122807 735.8141    5086
Urbana  CS  1   0.972332    4694.035088 395.2647    4639
Urbana  CS  2   0.952922    5901.052632 570.4148    4781
Urbana  CS  3   0.804431    4303.859649 458.0383    4934
Urbana  CS  4   0.742634    6768.421053 535.3851    4857
Urbana  SCS 1   1.195837    6159.298246 678.0293    4537
Urbana  SCS 2   1.267285    6090.526316 670.7419    4890
Urbana  SCS 3   1.08571 4939.649123 335.2923    4842
Urbana  SCS 4   1.20097 5262.45614  562.5674    4608
")

I would like to see if any of the variables (including factors in my treatment have significant effect on yield. So here, I am running stepwise regression which determined the best fit glm(formula = yield ~ treatment + Location, family = "gaussian", data = df). However, I am not sure why it is not showing the results from summary(glm.step) for CC and Urbana from variables treatment and Location. Can someone please explain me the reason?

n = nrow(df)
null = glm(yield ~ 1, data = df, family = "gaussian")
full = glm(
  yield ~ .,
  family = "gaussian",
  data = df
)

glm.step <-
  step(
    null,
    scope = list(lower = null, upper = full),
    direction = "forward",
    k = log(n),
    trace = 0
  )
glm.step
summary(glm.step)

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ Hopefully, the current answer by @Patrick provides an answer to your question. But I'll add that the results you may be seeking might be produced by e.g. library(car); Anova(Model), where Model is your final glm model. From there, you can use the emmeans package to further compare the treatments or interactions of treatments. $\endgroup$ Oct 19, 2019 at 16:04

1 Answer 1

3
$\begingroup$

R is using "LocationMonmouth" and "TreatmentCC" as reference levels and the other levels are given as relative to that.

By way of example let's consider a simple model that only has the location variable in it:

What you want to see is:

$yield = intercept + Monmouth*IndicatorMonmouth + Urbana*IndicatorUrbana$

Where $IndicatorMonmouth$ is 1 if the data point is from Monmouth and 0 if it's not.

But as you can see we have 3 variables to solve for now and our data only has two levels in it. So we have more unknowns than knowns and thus it is impossible to solve that formulation of the model. So what R and glm() is doing for you is rephrasing your model as follows:

$yield = intercept + (Urbana - Monmouth)*IndicatorUrbana$

Now there are only two parameters to estimate: intercept and (Urbana - Monmouth). So it has chosen to put the Monmouth data into the intercept term and present you with a parameter that is actually the difference between Urbana and Monmouth. Unfortunately R doesn't rename the variables or explicitly tell you it did that.

Same thing applies for TreatmentCC. I hope this also clears up how to interpret the output from the glm function.

Reference: https://stats.idre.ucla.edu/r/modules/coding-for-categorical-variables-in-regression-models/

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.