تمرین درس R(رگرسیون غیر خطی در R) هلیا علیزده

هلیا علیزاده تمرین درس R رگرسیون غیر خطی در R  

• برای بررسی معادله روابط متغییر هایی که نوه ی تغییرات متقابل انها به صورت خطی نباشند از رگرسیون غیر خطی استفاده می شود .

رگرسیون غیر خطی انواع خاصی دارد که رد ان یک مدل ریاضیتی غیر خطی برای تشریح رابطه متغییر ها ب کار گرفته می شود . در مدل غیر خطی  حداقل یکی از ضرایب مدل در شکل غیر خطی حضور دارد. یکی از مشهورترین مدل های غیر خطی مدل های رشد نمایی است . تجزیه رگرسیون غیر خطی در نرم افزار R  با تابع nls() صورت میگیرد . شکل کلی استفاده از تابع  غیرخطی به صورت زیر است : ((رابطه غیر خطی)~  متغیر وابسته ) nls نحوه ی استفاده از تابع nls مثال : رابطه وزن جوجه ها بر حسب کیلو گرم (آبجکت weight)  با سن آنها بر حسب روز (آبجکت time) با اندازه گیری یازده  فرد مورد مطالعه قرار گرفته است . برای یافتن معادلات غیر خطی و رسم منحنی آنها به مقادیر اولیه ضرایب (پارامتر ها) مدل نیاز است . ابتدا نمودار پراکنش نقاط ترسیم می شود تا رابطه دو متغییر بیان گردد . در همین راستا نیز می توانید از مقاله ای تخصصی رگرسیون ریچ را بررسی کنید.

> time<-c(6,7,8,9,10,11,12,13,14,15,16)

> weight<-c(29,52,79,125,181,261,425,738,1130,1882,2812)

> plot(time,weight)

  در نمودار مشاهده می کنیم که تغییرات وزن جوجه ها در هر زمان ب صورت یک منحنی (غیر خطی) افزایش یافته است . استفاده از رگرسیمو غیر خطی دقت بیشتری در تخمین وزن جوجه دارد   استفاده از تابع nls

> para.st<-c(a=0.006,b=4.600)

> nlmodel1<-nls(weight~(a*time^b),start=para.st,trace=T)

693847.5 :  0.006 4.600

683618.5 :  0.005566805 4.627912047

673420.9 :  0.005169855 4.655487761

663254.8 :  0.004805834 4.682726097

653120.7 :  0.004471749 4.709626300

643019.7 :  0.004164896 4.736187896

632953.4 :  0.003882832 4.762410665

632450 :  0.00336386 4.81417864

630663.7 :  0.002924614 4.864847248

627549.1 :  0.002551942 4.914382422

623094.7 :  0.002234944 4.962756063

617317.7 :  0.001964585 5.009945907

610259.6 :  0.001733368 5.055935377

601981.8 :  0.001535064 5.100713305

592561.8 :  0.001364497 5.144273649

582088.8 :  0.001217356 5.186615153

570660.7 :  0.001090046 5.227740969

558380.8 :  0.0009795636 5.2676582659

545355.1 :  0.0008833981 5.3063778449

531690 :  0.000799443 5.343913729

517490.4 :  0.0007259293 5.3802827905

502858.2 :  0.0006613679 5.4155043504

487890.7 :  0.0006045025 5.4495998273

472680.3 :  0.000554271 5.482592408

471126.9 :  0.000465275 5.546421017

462810.3 :  0.0003953606 5.6064939036

448930.8 :  0.0003399161 5.6628571207

430723.7 :  0.0002955267 5.7156051004

409361.6 :  0.0002596515 5.7648688326

385896 :  0.0002303887 5.8108050259

361228 :  0.0002063054 5.8535867796

351469.4 :  0.0001663235 5.9332050763

322507.9 :  0.0001387286 6.0024949047

284170.7 :  0.0001191369 6.0623167246

261374.8 :  9.055774e-05 6.165176e+00

241809.6 :  0.0000603775 6.3158220281

20093.78 :  4.643976e-05 6.454989e+00

13958.51 :  4.889906e-05 6.443546e+00

13941.47 :  4.898083e-05 6.443344e+00

13941.47 :  4.898192e-05 6.443336e+00

> summary(nlmodel1)

Formula: weight ~ (a * time^b)

Parameters:

   Estimate Std. Error t value Pr(>|t|)   

a 4.898e-05  2.043e-05   2.397   0.0401 * 

b 6.443e+00  1.529e-01  42.137 1.19e-11 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 39.36 on 9 degrees of freedom

Number of iterations to convergence: 39

Achieved convergence tolerance: 6.912e-07

  همانگونه که می شود از براطش مدل  رابطه وزن جوجه و سن آن بصورت    به دست آمده است . برای رسم نمودار پراکنش نقاط و مشاهده منحنی دستور زیر را به کار می بریم

(plot(time,weight<

((lines(time,fitted.values(nlmodel1<

  مقادیر برآورد اولیه ضرایب مدل (b,a) ب صورت تجربی و یا با ملاحضات تئوری توصیف می شود در این مثال از براورد های اولیه a=0.006  و b=4.6 استفاده گردید . روش ساده برآورد مقادیر اولیه ضرایب مدل استفاده از رگرسیون خطی یا داده های تبدیل شده می باشد . در مثال زیر نحوه برآورد این مقادیر اولیه در مورد این داده ها نشان داده شده است .

(( lm(log(weight)~log(time<

 

Call:

((lm(formula = log(weight) ~ log(time

Coefficients:

(Intercept)    log(time) 

     -5.174        4.611 

  همانگونه که مشاهده می شود مدل تبدیل شده بصورت زیر می باشد

(ln(weight)=-5.174+4.611 ln(time

b=4.6 و  برآورد های اولیه ضرایب مدل نمایی می باشد . در مواردی که مدل ها ناشناخته اند و برآورد اولیه پارامترها کمی دشوار باشد می توان از تابع getInitial() و یک تابع همراه دیگر (بر حسب نوع رابطه ) استفاده کرد . در مثال زیرداده هایی که دارای رابطه به شکل منحنی لجستیک هستند مورد استفاده قرارگرفته است تابه کارگیری توابع  nls() و getInitial() و Sslogis()  مدل غیر خطی مناسب برازش داده شود .   داده های پیشرفت شدت بیماری (severity) در طی زمان (time) بصورت داده چارچوب دار (با استفاده از تابع data.frame()) در ابجکت disease نگهداری شده اند . ایتدا نمودار پراکنش نقاط با استفاده از تابع plot() ترسیم می شود  

 disease<

data.frame(severity=c(42,51,59,64,76,93,106,125,149,171,199),time=c(0,1,2,3,4,5,6,7,8,9,10))

( plot(disease$time,disease$severity<

  همانگونه که مشاهده می شود نوع ارتباط 2متغیر از نوع غیر خطی است . تجربیات قبلی  و شکل پرا کنش نقاط نشان می دهد که نوع رابطه بصورت یک تابع لجستیک  می باشد . ابتدا با استفاده از توابع  getInitial() و Sslogis() مقادیر اولیه ضرایب مدل لجستیک تخمین زده می شود .

> getInitial(severity~SSlogis(time,a,xmid,scale),data=disease)

          a        xmid       scale

9907.723356   34.581254    6.329369

تخمین اولیه ضریب b نیز برابر با کسر xmid/scale است . تخمین های اولیه را درون آبجکت para1  نگهداری می نماییم .

> para1<-c(a=9907.73,b=34.58/6.33,r=1.633)

> para1

          a           b           r

9907.730000    5.462875    1.633000

  مقادیر تخمین های اولیه را رد برهان start= قرار می دهیم .  

> nlmodel2<-nls(severity~(a/(1+exp(b-r*time))),data=disease,start=para1)

> nlmodel2

Nonlinear regression model

  model: severity ~ (a/(1 + exp(b - r * time)))

   data: disease

       a        b        r

9902.828    5.463    0.158

 residual sum-of-squares: 32.39

 

Number of iterations to convergence: 31

Achieved convergence tolerance: 2.455e-06

  چنانچه بخواهیم دوره های تخمین ضرایب را مشاهده کنیم برهان trace=T را نیز در تابع فوق قرار می دهیم . همانگونه که در خروجی فوق مشاهده می شود رابطه شدت بیماری (severity) و زمان (time) بصورت  می باشد .  

> summary(nlmodel2)

Formula: severity ~ (a/(1 + exp(b - r * time)))

Parameters:

   Estimate Std. Error t value Pr(>|t|)   

a 9.903e+03  5.499e+04   0.180    0.862   

b 5.463e+00  5.556e+00   0.983    0.354   

r 1.580e-01  1.082e-02  14.606 4.73e-07 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.012 on 8 degrees of freedom

Number of iterations to convergence: 31

Achieved convergence tolerance: 2.455e-06

    برای مشاهده منحنی این رابطه به همراه نقاط مشاهده شده دستورات زیر را اضافه می نماییم .  

> plot(disease$time,disease$severity)

> curve(9902.134/(1+exp(5.463-0.158*x)),from=0,to=10,add=TRUE)

  همانگونه که مشاهده می شود با استفاده از  تابع curve() و معادله منحنی ، شکل آن ترسیم شده است برهان add=TRUE به شما اجازه می دهد تا منحنی را روی نمودار پراکنش نقاط اضافه نمایید . در این مورد که رابطه متغییر ها ب صورت تابع لجستیک بود از تابع Sslogis() برای تخمین اولیه ضرایب مدل استفاده گردید . در نرم افزار R توابع مشابهی برای تخمین ضرایب مدل های نمایی همانند Ssexp()  و مدل ویبولی  همانند Ssweibull() وجود دارد .   در زیر مثالی برای تابع nls() گفته ایم :

Examples

require(graphics)

DNase1 <- subset(DNase, Run == 1)

 

## using a selfStart model## با استفاده از مدل خود شروع

fm1DNase1 <- nls(density ~ SSlogis(log(conc), Asym, xmid, scal), DNase1)

summary(fm1DNase1)

## the coefficients only:  # # فقط ضرایب:

coef(fm1DNase1)

## including their SE, etc:  ## جمله SE خود، و غیره:

coef(summary(fm1DNase1))

 

## using conditional linearity  ## با استفاده از خطی مشروط

fm2DNase1 <- nls(density ~ 1/(1 + exp((xmid - log(conc))/scal)),

                 data = DNase1,

                 start = list(xmid = 0, scal = 1),

                 algorithm = "plinear")

summary(fm2DNase1)

 

## without conditional linearity  ## با استفاده از خطی مشروط

fm3DNase1 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)),

                 data = DNase1,

                 start = list(Asym = 3, xmid = 0, scal = 1))

summary(fm3DNase1)

## using Port's nl2sol algorithm   ## با استفاده از الگوریتم nl2sol پورتر

fm4DNase1 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)),

                 data = DNase1,

                 start = list(Asym = 3, xmid = 0, scal = 1),

                 algorithm = "port")

summary(fm4DNase1)

## weighted nonlinear regression  ## رگرسیون غیر خطی وزن دار

Treated <- Puromycin[Puromycin$state == "treated", ]

weighted.MM <- function(resp, conc, Vm, K)

{

    ## Purpose: exactly as white book p. 451 -- RHS for nls()

    ##  Weighted version of Michaelis-Menten model

    ## ----------------------------------------------------------

    ## Arguments: 'y', 'x' and the two parameters (see book)

    ## ----------------------------------------------------------

    ## Author: Martin Maechler, Date: 23 Mar 2001

## هدف: دقیقا به عنوان سفید ص کتاب. 451 - RHS برای NLS () ## نسخه وزنی مدل Michaelis به-منتن ## ------------------------------------------------ ---------- برهان ##: 'Y'، 'X' و دو پارامتر (کتاب را ببینید) ## ------------------------------------------------ ---------- ## نویسنده: مارتین Maechler، تاریخ: 2001 مارس 23

    pred <- (Vm * conc)/(K + conc)

    (resp - pred) / sqrt(pred)

}   Pur.wt <- nls( ~ weighted.MM(rate, conc, Vm, K), data = Treated, start = list(Vm = 200, K = 0.1))

summary(Pur.wt)

## Passing arguments using a list that can not be coerced to a data.frame

## عبور استدلال با استفاده از یک لیست است که می تواند به یک data.frame تومی نمیتواند مجبور

lisTreat <- with(Treated,

                 list(conc1 = conc[1], conc.1 = conc[-1], rate = rate))

weighted.MM1 <- function(resp, conc1, conc.1, Vm, K)

{

     conc <- c(conc1, conc.1)

     pred <- (Vm * conc)/(K + conc)

    (resp - pred) / sqrt(pred)

}

Pur.wt1 <- nls( ~ weighted.MM1(rate, conc1, conc.1, Vm, K),

               data = lisTreat, start = list(Vm = 200, K = 0.1))

stopifnot(all.equal(coef(Pur.wt), coef(Pur.wt1)))

## Chambers and Hastie (1992) Statistical Models in S  (p. 537):

## If the value of the right side [of formula] has an attribute called

## 'gradient' this should be a matrix with the number of rows equal

## to the length of the response and one column for each parameter.

## اتاق و HASTIE (1992) مدل های آماری در S (ص 537): ## اگر مقدار سمت راست [از فرمول] دارای یک ویژگی به نام ##، شیب، این باید یک ماتریس با تعداد ردیف برابر باشد # # به طول پاسخ و یک ستون برای هر پارامتر.

weighted.MM.grad <- function(resp, conc1, conc.1, Vm, K)

{

  conc <- c(conc1, conc.1)

  K.conc <- K+conc

  dy.dV <- conc/K.conc

  dy.dK <- -Vm*dy.dV/K.conc

  pred <- Vm*dy.dV

  pred.5 <- sqrt(pred)

  dev <- (resp - pred) / pred.5

  Ddev <- -0.5*(resp+pred)/(pred.5*pred)

  attr(dev, "gradient") <- Ddev * cbind(Vm = dy.dV, K = dy.dK)

  dev

}

Pur.wt.grad <- nls( ~ weighted.MM.grad(rate, conc1, conc.1, Vm, K),

                   data = lisTreat, start = list(Vm = 200, K = 0.1))

rbind(coef(Pur.wt), coef(Pur.wt1), coef(Pur.wt.grad))

## In this example, there seems no advantage to providing the gradient.

## In other cases, there might be.

## The two examples below show that you can fit a model to

## artificial data with noise but not to artificial data

## without noise.

## در این مثال، به نظر می رسد هیچ مزیت به ارائه شیب وجود دارد. ## در موارد دیگر، ممکن است وجود داشته باشد. # # دو مثال زیر نشان می دهد که شما می توانید یک مدل را به تناسب ## داده های مصنوعی با سر و صدا اما نه به داده های مصنوعی ## بدون سر و صدا.

x <- 1:10

y <- 2*x + 3                            # perfect fit

yeps <- y + rnorm(length(y), sd = 0.01) # added noise

nls(yeps ~ a + b*x, start = list(a = 0.12345, b = 0.54321))

## terminates in an error, because convergence cannot be confirmed:

## خاتمه یافته در یک خطا، به دلیل همگرایی نمی توان تایید کرد:

try(nls(y ~ a + b*x, start = list(a = 0.12345, b = 0.54321)))

## the nls() internal cheap guess for starting values can be sufficient:

# # NLS به () حدس ارزان داخلی برای شروع مقادیر می تواند به اندازه کافی:

x <- -(1:100)/10

y <- 100 + 10 * exp(x / 2) + rnorm(x)/10

nlmod <- nls(y ~  Const + A * exp(B * x))

plot(x,y, main = "nls(*), data, true function and fit, n=100")

curve(100 + 10 * exp(x / 2), col = 4, add = TRUE)

lines(x, predict(nlmod), col = 2)

## The muscle dataset in MASS is from an experiment on muscle

## contraction on 21 animals.  The observed variables are Strip

## (identifier of muscle), Conc (Cacl concentration) and Length

## (resulting length of muscle section).

## مجموعه داده عضلانی در جرم از یک آزمایش در مورد عضله ## انقباض در 21 حیوانات. متغیرهای مشاهده شده است نوار ## (شناسه عضلانی)، بتنی (غلظت CACL) و طول ## (طول بخش عضله ناشی).

utils::data(muscle, package = "MASS")

## The non linear model considered is

##       Length = alpha + beta*exp(-Conc/theta) + error

## where theta is constant but alpha and beta may vary with Strip.

## مدل غیر خطی در نظر گرفته است ## طول = آلفا + بتا * EXP (-Conc / تتا) + خطا ## که در آن تتا ثابت است اما آلفا و بتا ممکن است با نوار متفاوت است.

with(muscle, table(Strip)) # 2, 3 or 4 obs per strip

## We first use the plinear algorithm to fit an overall model,

## ignoring that alpha and beta might vary with Strip.

## ما برای اولین بار استفاده از الگوریتم خطی به جا یک مدل به طور کلی، ## نادیده گرفتن که آلفا و بتا ممکن است با نوار متفاوت است.

musc.1 <- nls(Length ~ cbind(1, exp(-Conc/th)), muscle,

              start = list(th = 1), algorithm = "plinear")

summary(musc.1)

## Then we use nls' indexing feature for parameters in non-linear

## models to use the conventional algorithm to fit a model in which

## alpha and beta vary with Strip.  The starting values are provided

## by the previously fitted model.

## Note that with indexed parameters, the starting values must be

## given in a list (with names):

## سپس با استفاده از ویژگی های نمایه سازی NLS برای پارامترها در غیر خطی مدل ## به استفاده از الگوریتم معمولی به تناسب مدل که در آن ## آلفا و بتا با نوار متفاوت است. ارزش شروع ارائه شده است ## با استفاده از مدل قبلا نصب شده است. ## توجه داشته باشید که با پارامترهای نمایه، ارزش شروع باید ## داده شده در یک لیست (با نام):

b <- coef(musc.1)

musc.2 <- nls(Length ~ a[Strip] + b[Strip]*exp(-Conc/th), muscle,

              start = list(a = rep(b[2], 21), b = rep(b[3], 21), th = b[1]))

summary(musc.2)