library(equatiomatic)
library(readxl)

Question 1

Part (a)

spider_data <- read_excel("spider_data.xlsx")
#pdead = spider_data$pdead
#dose = spider_data$dose
#male = (spider_data$sex)
#alive = spider_data$alive
model = lm(formula = pdead ~ dose + male, data = spider_data);
summary(model)

Call:
lm(formula = pdead ~ dose + male, data = spider_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.15043 -0.07913 -0.01601  0.07127  0.17125 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.04578    0.09130  -0.501    0.631    
dose         0.75225    0.08464   8.888 4.63e-05 ***
male         0.06783    0.08389   0.809    0.445    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1326 on 7 degrees of freedom
Multiple R-squared:  0.9192,    Adjusted R-squared:  0.8961 
F-statistic: 39.82 on 2 and 7 DF,  p-value: 0.0001499

Use the equatiomatic package to display the model with or without coefficients!

extract_eq(model,use_coefs = FALSE) 

\[ \operatorname{pdead} = \alpha + \beta_{1}(\operatorname{dose}) + \beta_{2}(\operatorname{male}) + \epsilon \]

extract_eq(model,use_coefs = TRUE) 

\[ \operatorname{pdead} = -0.05 + 0.75(\operatorname{dose}) + 0.07(\operatorname{male}) + \epsilon \]

We can easily work out the answer to part a) using these:

## At 0.9 mmol
(-0.05 + 0.9 * 0.75 + 0.07 * 1) * 100 # male
[1] 69.5
(-0.05 + 0.9 * 0.75 + 0.07 * 0) * 100 # female
[1] 62.5
## At 1.4 mmol
(-0.05 + 1.4 * 0.75 + 0.07 * 1) * 100 # male
[1] 107
(-0.05 + 1.4 * 0.75 + 0.07 * 0) * 100 # female
[1] 100

Note how some of these values are impossible!

Part (b)

model2 <- glm(formula = pdead ~ dose + male, family = quasibinomial(link="logit"), data = spider_data);
summary(model2)

Call:
glm(formula = pdead ~ dose + male, family = quasibinomial(link = "logit"), 
    data = spider_data)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.18266  -0.04974  -0.01073   0.05891   0.12985  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -4.7322     0.3678 -12.866 3.98e-06 ***
dose          6.7843     0.4761  14.249 1.99e-06 ***
male          0.7744     0.2309   3.354   0.0122 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.01119487)

    Null deviance: 7.918937  on 9  degrees of freedom
Residual deviance: 0.075013  on 7  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 7

The model takes 10 datapoints for each sex. The null model has one parameter, and so has 9 DoF. The fitted model has 3 parameters, and so has 7 DoF.

Part (c)

extract_eq(model2, use_coefs = TRUE) 

\[ \log\left[ \frac { P( \operatorname{pdead} = \operatorname{1} ) }{ 1 - P( \operatorname{pdead} = \operatorname{1} ) } \right] = -4.73 + 6.78(\operatorname{dose}) + 0.77(\operatorname{male}) + \epsilon \]

You can similarly find the linear predictors by plugging the numbers from the question sheet into the formula above (note I don’t know why the question sheet answers are slightly different, I suspect he has used a larger dataset than than what is given.)

Part (d)

It just so happens that for a dose of 0.7021, the linear estimator is close to zero. Plugging this into the inverse logit gives a value of 0.5 for the proportion of dead female spiders at this concentration. Similarly for male spiders, one finds that we expect the proportion of dead spiders at this concentration to be around 0.7. Taking the ratio of these two values shows that male spiders are 1.4 times more likely to die than female spiders at this concentration.

Question 2

Showing that the maximum liklihood estimate of the probability in a binomial distribution if one gets \(y\) successes in \(N\) trials is given by \(p=y/n\):

The PMF of the binomial distribution is given as: \[B(n,p) = {n \choose k} p^k (1-p)^{n-k}\] Taking logarithms:

\[\ln B(n,p) = \ln {n \choose k} + k \ln p + (n-k) \ln (1-p) \]

We would like to find the probability which maximises the liklihood, this is the same as taking \(\frac{d\ln B(n,p)}{dp} = 0\):

\[\frac{d\ln B(n,p)}{dp} = \frac{k}{p} - \frac{n-k}{1-p} = 0\] \[ \frac{k}{p} = \frac{n-k}{1-p} \rightarrow \frac{1-p}{p} = \frac{n-k}{k} \rightarrow \frac{1}{p} = \frac{n}{k}\]

Therefore we can see that:

\[p = \frac{k}{n}\]

Question 3 and 4

Unfortunately the raw data for these questions are not available! See the worked solutions uploaded by the lecturer!

LS0tCnRpdGxlOiAiR0xNIHN1cGVydmlzaW9uIgphdXRob3I6ICJHZW9yZ2VvcyBIYXJkbyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKCgpgYGB7cn0KbGlicmFyeShlcXVhdGlvbWF0aWMpCmxpYnJhcnkocmVhZHhsKQpgYGAKCiMjIFF1ZXN0aW9uIDEKCiMjIyBQYXJ0IChhKQoKYGBge3J9CnNwaWRlcl9kYXRhIDwtIHJlYWRfZXhjZWwoInNwaWRlcl9kYXRhLnhsc3giKQpgYGAKCmBgYHtyfQojcGRlYWQgPSBzcGlkZXJfZGF0YSRwZGVhZAojZG9zZSA9IHNwaWRlcl9kYXRhJGRvc2UKI21hbGUgPSAoc3BpZGVyX2RhdGEkc2V4KQojYWxpdmUgPSBzcGlkZXJfZGF0YSRhbGl2ZQptb2RlbCA9IGxtKGZvcm11bGEgPSBwZGVhZCB+IGRvc2UgKyBtYWxlLCBkYXRhID0gc3BpZGVyX2RhdGEpOwpzdW1tYXJ5KG1vZGVsKQpgYGAKVXNlIHRoZSBlcXVhdGlvbWF0aWMgcGFja2FnZSB0byBkaXNwbGF5IHRoZSBtb2RlbCB3aXRoIG9yIHdpdGhvdXQgY29lZmZpY2llbnRzIQpgYGB7ciwgcmVzdWx0cyA9ICJhc2lzIn0KZXh0cmFjdF9lcShtb2RlbCx1c2VfY29lZnMgPSBGQUxTRSkgCmV4dHJhY3RfZXEobW9kZWwsdXNlX2NvZWZzID0gVFJVRSkgCmBgYApXZSBjYW4gZWFzaWx5IHdvcmsgb3V0IHRoZSBhbnN3ZXIgdG8gcGFydCBhKSB1c2luZyB0aGVzZToKCmBgYHtyfQojIyBBdCAwLjkgbW1vbAooLTAuMDUgKyAwLjkgKiAwLjc1ICsgMC4wNyAqIDEpICogMTAwICMgbWFsZQooLTAuMDUgKyAwLjkgKiAwLjc1ICsgMC4wNyAqIDApICogMTAwICMgZmVtYWxlCiMjIEF0IDEuNCBtbW9sCigtMC4wNSArIDEuNCAqIDAuNzUgKyAwLjA3ICogMSkgKiAxMDAgIyBtYWxlCigtMC4wNSArIDEuNCAqIDAuNzUgKyAwLjA3ICogMCkgKiAxMDAgIyBmZW1hbGUKYGBgCk5vdGUgaG93IHNvbWUgb2YgdGhlc2UgdmFsdWVzIGFyZSBpbXBvc3NpYmxlIQoKIyMjIFBhcnQgKGIpCgpgYGB7cn0KbW9kZWwyIDwtIGdsbShmb3JtdWxhID0gcGRlYWQgfiBkb3NlICsgbWFsZSwgZmFtaWx5ID0gcXVhc2liaW5vbWlhbChsaW5rPSJsb2dpdCIpLCBkYXRhID0gc3BpZGVyX2RhdGEpOwpzdW1tYXJ5KG1vZGVsMikKYGBgCgpUaGUgbW9kZWwgdGFrZXMgMTAgZGF0YXBvaW50cyBmb3IgZWFjaCBzZXguIFRoZSBudWxsIG1vZGVsIGhhcyBvbmUgcGFyYW1ldGVyLCBhbmQgc28gaGFzIDkgRG9GLiBUaGUgZml0dGVkIG1vZGVsIGhhcyAzIHBhcmFtZXRlcnMsIGFuZCBzbyBoYXMgNyBEb0YuCgojIyMgUGFydCAoYykKCmBgYHtyLCByZXN1bHRzID0gImFzaXMifQpleHRyYWN0X2VxKG1vZGVsMiwgdXNlX2NvZWZzID0gVFJVRSkgCmBgYApZb3UgY2FuIHNpbWlsYXJseSBmaW5kIHRoZSBsaW5lYXIgcHJlZGljdG9ycyBieSBwbHVnZ2luZyB0aGUgbnVtYmVycyBmcm9tIHRoZSBxdWVzdGlvbiBzaGVldCBpbnRvIHRoZSBmb3JtdWxhIGFib3ZlIChub3RlIEkgZG9uJ3Qga25vdyB3aHkgdGhlIHF1ZXN0aW9uIHNoZWV0IGFuc3dlcnMgYXJlIHNsaWdodGx5IGRpZmZlcmVudCwgSSBzdXNwZWN0IGhlIGhhcyB1c2VkIGEgbGFyZ2VyIGRhdGFzZXQgdGhhbiB0aGFuIHdoYXQgaXMgZ2l2ZW4uKQoKIyMjIFBhcnQgKGQpCgpJdCBqdXN0IHNvIGhhcHBlbnMgdGhhdCBmb3IgYSBkb3NlIG9mIDAuNzAyMSwgdGhlIGxpbmVhciBlc3RpbWF0b3IgaXMgY2xvc2UgdG8gemVyby4gUGx1Z2dpbmcgdGhpcyBpbnRvIHRoZSBpbnZlcnNlIGxvZ2l0IGdpdmVzIGEgdmFsdWUgb2YgMC41IGZvciB0aGUgcHJvcG9ydGlvbiBvZiBkZWFkIGZlbWFsZSBzcGlkZXJzIGF0IHRoaXMgY29uY2VudHJhdGlvbi4gU2ltaWxhcmx5IGZvciBtYWxlIHNwaWRlcnMsIG9uZSBmaW5kcyB0aGF0IHdlIGV4cGVjdCB0aGUgcHJvcG9ydGlvbiBvZiBkZWFkIHNwaWRlcnMgYXQgdGhpcyBjb25jZW50cmF0aW9uIHRvIGJlIGFyb3VuZCAwLjcuIFRha2luZyB0aGUgcmF0aW8gb2YgdGhlc2UgdHdvIHZhbHVlcyBzaG93cyB0aGF0IG1hbGUgc3BpZGVycyBhcmUgMS40IHRpbWVzIG1vcmUgbGlrZWx5IHRvIGRpZSB0aGFuIGZlbWFsZSBzcGlkZXJzIGF0IHRoaXMgY29uY2VudHJhdGlvbi4gCgojIyBRdWVzdGlvbiAyCgpTaG93aW5nIHRoYXQgdGhlIG1heGltdW0gbGlrbGlob29kIGVzdGltYXRlIG9mIHRoZSBwcm9iYWJpbGl0eSBpbiBhIGJpbm9taWFsIGRpc3RyaWJ1dGlvbiBpZiBvbmUgZ2V0cyAkeSQgc3VjY2Vzc2VzIGluICROJCB0cmlhbHMgaXMgZ2l2ZW4gYnkgJHA9eS9uJDoKClRoZSBQTUYgb2YgdGhlIGJpbm9taWFsIGRpc3RyaWJ1dGlvbiBpcyBnaXZlbiBhczoKJCRCKG4scCkgPSB7biBcY2hvb3NlIGt9IHBeayAoMS1wKV57bi1rfSQkClRha2luZyBsb2dhcml0aG1zOiAKCiQkXGxuIEIobixwKSA9IFxsbiB7biBcY2hvb3NlIGt9ICsgayBcbG4gcCArIChuLWspIFxsbiAoMS1wKSAkJAoKV2Ugd291bGQgbGlrZSB0byBmaW5kIHRoZSBwcm9iYWJpbGl0eSB3aGljaCBtYXhpbWlzZXMgdGhlIGxpa2xpaG9vZCwgdGhpcyBpcyB0aGUgc2FtZSBhcyB0YWtpbmcgJFxmcmFje2RcbG4gQihuLHApfXtkcH0gPSAwJDoKCiQkXGZyYWN7ZFxsbiBCKG4scCl9e2RwfSA9IFxmcmFje2t9e3B9IC0gXGZyYWN7bi1rfXsxLXB9ID0gMCQkCiQkIFxmcmFje2t9e3B9ID0gXGZyYWN7bi1rfXsxLXB9ICBccmlnaHRhcnJvdyBcZnJhY3sxLXB9e3B9ID0gXGZyYWN7bi1rfXtrfSBccmlnaHRhcnJvdyBcZnJhY3sxfXtwfSA9IFxmcmFje259e2t9JCQKClRoZXJlZm9yZSB3ZSBjYW4gc2VlIHRoYXQ6CgokJHAgPSBcZnJhY3trfXtufSQkCgojIyBRdWVzdGlvbiAzIGFuZCA0CgpVbmZvcnR1bmF0ZWx5IHRoZSByYXcgZGF0YSBmb3IgdGhlc2UgcXVlc3Rpb25zIGFyZSBub3QgYXZhaWxhYmxlISBTZWUgdGhlIHdvcmtlZCBzb2x1dGlvbnMgdXBsb2FkZWQgYnkgdGhlIGxlY3R1cmVyIQ==