R language GLM generalized linear model: logistic regression, Poisson regression fitting mouse clinical trial data (dose and response) examples and self-test questions

Link to the original text ：http://tecdat.cn/?p=27464 

The source of the original text is ： The official account of the tribal public

 

In fitting GLM（ And check the residuals ） after , have access to z Test - test the significance of the estimated parameters , Compare the estimated value with its standard error .

GLM Examples of model fitting and analysis

Example 1. Mouse data GLM modeling （ Dose and response ）

a) We input data and fit logistic regression models .

 

> summary(it1.lt)

 

1-pchisq(17.6,24) 

 

Model ：

It can be compared with the complete model . And deviation value 17.639 dependent P value 0.82（> 0.10） This means that there is no significant evidence to reject the fitted model .

anova(fi.lgi)

 

 

1-pchisq(35.8-17.69, 25-24)

 

The empty model does not contain predictive variables , stay 25 A degree of freedom (df) The deviation on is 35.89. When the covariate x When adding to an empty model , The variation of deviation is 35.890-17.639=18.25. And the degree of freedom is 25-24=1 Chi square distribution of , Its P The value is 1.93 × 10 -5 Very remarkable .

Therefore, the model cannot be deleted x To simplify the .x Of the coefficient of t Inspection is also important （P value 0.0065<0.01）.

Intercept ？ Can you delete it ？

> plotx, itte(fi1log,typ"
> pot(,y

chart 1： Logistic regression of mouse data and fitting values .
b） We fit a model with probabilistic links .

> summary()

 

Supporting model ：

Again , Both of these parameters are important （P value <0.01）

> anova

 

> 1 - pchisq(35.89-17.49 25-24)

> lines(x, fitte

  add to x The variation of time deviation is significant （P value =）.
The model cannot be deleted x To simplify the .

chart 2： Mouse data and fitting values （ Dotted line ： Probabilistic Links ）.
Use probit Linked models are slightly better than using logit Linked models , Because the deviation is smaller . In two models ,x The coefficients of are significant （P value <0.01）, This means that the effect increases as the dose increases .
 

Example 2. Clinical trial data （ Dose and response ） Of GLM modeling .

a) We input data , Then fit the logistic regression model
 

> summary(it2.it)

 

 1-pchisq(13.63,6)

 

And deviation value 13.633 dependent P The value is 0.034<0.05.5% The level of rejection fits the model .
in the light of x Plotting residuals reveals a dependency pattern .

> plot(x, reid(it2.it))

 

chart 3： Only with x The residual diagram of the fitting model .

plt(fitdft2lit reid(fi2lot))

  So we're going to x2 Add to model .

> summary(ft2qlt)

 

> 1-pchisq(5.1, 5)

Deviation from 13.633 Reduced to 5.107, Not significant （P value =0.403>0.05）.
therefore , We cannot reject this model through evidence of bias .

plot(fitted(fit2logit), resid(.logit))

 

chart 4： with x2 The residual diagram of the fitting model .

The residuals now look random .
The fitting model is

And all parameter estimates are significant （5%）. Logarithmic probability Depend on in a quadratic way x.
 

  Example 3. AIDS data , Poisson

a) We input the data and use the default logarithmic link to fit the Poisson regression model .

> smary(fit.lg)

 

> 1-pchisq

 

deviation 29.654 Of P The value is 0.005<0.01 ⇒ Model rejected .

plot(fit3esuals)

 

chart 5： Fitting model 3

Residual diagram of b) For year index x The residual graph of shows the dependency pattern . So we add .
 

> summary(fi3.lg)

> 1-pchisq
[1] 0.1279

16.371 The deviation of （P The value is 0.1279>0.10） Not significant . Fitting model

You can't refuse on the ground of deviation . But the residual plot shows only slightly more random patterns than before .
 

chart 6： Fitting model 3 And x2 Residual diagram of

The model can be improved by using non canonical links .

> summary(ft3st)

 

> 1 - pchisq(16.9, 12)
[1] 0.153

  Residual error of fitting model y = (-0.27571 +0.49277x)2 + e Shows a more random pattern .12 df The deviation in 16.905 Slightly higher than the previous model 16.371（df=11）, But still not significant （P value =0.1532>0.10）.AIC smaller , by 73.833<75.298. therefore , Models with square root links are preferred .

Constant items can be deleted （“ intercept ”） Do you ？

What other linking functions are available ？

Self test questions ：

Twenty tobacco budworm moths of each sex were exposed to different doses of the insecticide trans-cypermethrin. The numbers of budworm moths killed during a 3-day exposure were as follows for each sex (male, female) and dose level in mg’s.


Type the data into R as follows. Press Enter at the end of each line including blank lines.

num.killed <- scan()
1 4 9 13 18 20 0 2 6 10 12 16
sex <- scan()
0 0 0 0 0 0 1 1 1 1 1 1
dose <- scan()
1 2 4 8 16 32 1 2 4 8 16 32

Fit two models by doing the following.
 

ldose <- log(dose)/log(2) #convert to base-2 log dose
ldose #have a look
y <- cbind(num.killed, 20-num.killed) #add number survived
fit1 <- glm(y ~ ldose * sex, family=binomial(link=probit))
fit2 <- glm(y ~ sex + ldose, family=binomial(link=probit))

You may also run the following lines and refer to the chi-square distribution table

anova(fit1,test="Chisq")
summary(fit2)

1. What model is fitted in fit1? Write it formally and define all the terms. 

2. How is the model in fit2 differ from that in fit1? 

3. Does the model in fit1 fit the data adequately? Use deviance to answer this question.
4. Can the model in fit1 be simplified to the model in fit2? Use change in deviance to answer
this question. 
5. Can sex be removed from the model in fit2? Use change in deviance to answer this ques
tion. 
6. What are the maximum likelihood estimates of the parameters of the additive model? What
are their standard errors? Test the significance of each parameter using its estimate and
standard error. 
7. How does the probability of a kill change with log dose and sex of the budworm moth accord
ing to the additive model? 

(a) Derive the survival function S(t) of a lifetime T » E xp(‚). Find ¡logS(t) and comment on it.
(b) Calculate the Kaplan-Meier estimate for each group in the following.
Treatment Group:
6,6,6,6*,7,9*,10,10*,11*,13,16,17*,19*,20*,22,23,25*,32*,32*,34*,35
Control Group (no treatment):
1,1,2,2,3,4,5,5,8,8,8,8,11,11,12,15,17,22,23
Note that * indicates right censored data.
(c) Use the log rank test to compare the two groups of lifetimes.
 

---

The most popular insights

1.R Language diversity Logistic Logical regression The application case

2. Panel smooth transfer regression (PSTR) Analyze the case and realize

3.matlab Partial least squares regression in （PLSR） And principal component regression （PCR）

4.R Language Poisson Poisson Regression model analysis case

5.R Language mixing effect logistic regression Logistic Model analysis of lung cancer

6.r In language LASSO Return to ,Ridge Ridge return and Elastic Net Model implementation

7.R Logical regression of language 、Naive Bayes Bayes 、 Decision tree 、 Random forest algorithm predicts heart disease

8.python Using linear regression to predict stock prices

9.R Language uses logical regression 、 Decision tree and random forest are used to classify and predict credit data sets