Workshop in R and GLMs:

1
Workshop in R and GLMs 4

• Diane Srivastava
• University of British Columbia
• srivast_at_zoology.ubc.ca

2
Exercise

• Fit the binomial glm survival sizetreat
• 2. Fit the bionomial glm parasitism sizetreat
• 3. Predict what size has 50 parasitism in
  treatment 0

3
Predicting size for p0.5, treat0

• Output from logistic regression with logit link
  predicted loge (p/1-p) abx
• So when p0.5, solve log(1)abx

4
What is equation for treat 0? treat 1?
Coefficients Estimate Std. Error z
value Pr(gtz) (Intercept) -2.38462
0.16780 -14.211 lt2e-16 size
0.76264 0.04638 16.442 lt2e-16 treat
0.28754 0.23155 1.242 0.214
sizetreat -0.09477 0.06357 -1.491
0.136
5
Rlecture.csv
3.12
6
Model simplification

• Parsimonious/ Logical sequence (e.g. highest
  order interactions first)
• 2. Stepwise sequence
• 3. Bayesian comparison of candidate models (not
  covered)

7
ANCOVA Difference between categories.
Constant, doesnt depend on size
Depends on size
sizetreat sig
sizetreat ns
12
10
8
Logit parasitism
Logit parasitism
6
4
2
0
0
2
4
6
Plant size
Plant size
8
Deletion tests

• How to change your model quickly
• model2lt-update(model1,.-sizetreat)
• How to do a deletion test
• anova(reduced model, full model, test"Chi")
• Test for interaction in logit parasitism ANCOVA
• If not sig, remove and continue. If sig, STOP!
• 2. Test covariate If not sig, remove and
  continue. If sig, put back and continue
• 3. Test main effect

9
Code for parasitism analysis
gt dslt-read.table(file.choose(), sep",",
headerTRUE) ds gt attach(ds) gt
parlt-cbind(parasitism, 100-parasitism) par gt
m1lt-glm(parsizetreat, datads,
familybinomial) gt summary(m1) gt m2lt-update(m1,
.-sizetreat) gt summary(m2) gt anova(m2,m1,
test"Chi") gt m3lt-update(m2, .-size) gt
anova(m3,m2, test"Chi") gt m3lt-update(m2,
.-treat) gt anova(m3,m2, test"Chi")
10
Context (often) matters!
What is the p-value for treat in sizetreat? tr
eat? Stepwise regression step(model)
11
Jump height (how high ball can be raised off the
ground)
8
11
10
9
Feet off ground
Total SS 11.11
12
(No Transcript)
13
F1,13
14
Why do you think weight is correlated with jump
height?
15
An idea Perhaps if we took two people of
identical height, the lighter one might actually
jump higher? Excess weight may reduce ability to
jump high
16
lighter
heavier
17
Tall people can jump higher
Heavy people often tall (tall people often heavy)

Height
Jump

-
Weight
People light for their height can jump a bit more
18
Species.txt
Rothamsted Park Grass experiment started in 1856
19
Exercise (species.txt)

• dianelt-read.table(file.choose(), headerT)
  diane attach(diane)
• Univariate trends
• plot(SpeciesBiomass)
• plot(SpeciespH)
• Combined trends
• plot(SpeciesBiomass, type"n")
• points(SpeciespH"high"BiomasspH"high")
• points(SpeciespH"mid"BiomasspH"mid",
  pch16)
• points(SpeciespH"low"BiomasspH"low",
  pch0)

20
Exercise (species.txt)

• 1. With a normal distribution, fit pHBiomass
• check model dignostics
• test interaction for significance
• 2. With a poisson distribution, fit pH Biomass
• check model dignostics
• test interaction for significance

21
Moral of the story Make sure you KNOW what you
are modelling!
22
Exercise (species.txt)

• 1. Fit glm SpeciespH, familygaussian
• 2. Test if low and mid pH have the same effect
• this is a planned comparison

23
Further reading

• Statistics An Introduction using R
• (M.J. Crawley, Wiley publishers)
• Extending the linear model with R
• (JJ Faraway, Chapman Hall/CRC)

24
Code for Species analysis
gt m1lt-glm(SpeciespHBiomass, familygaussian,
datadiane) gt summary(m1) gt m2lt-update(m1,
.-pHBiomass) gt anova(m2,m1, test"Chi") gt
par(mfrowc(2,2)) plot(m1) gt m3lt-glm(SpeciespHB
iomass, familypoisson, datadiane) gt
m4lt-update(m3, .-pHBiomass) gt anova(m4,m3,
test"Chi") gt par(mfrowc(2,2))
plot(m3) gtPHlt-(pH!"high")0 gt m5lt-glm(SpeciespH,
familygaussian, datadiane) gt m6lt-update(m5,
.-pHPH) gt anova(m6,m5, test"Chi")