# Nelineární regrese # přednáška 16. prosince 2009 # # (původní data doc. Janků upravena, aby lépe přiléhala) # rm(list=ls()) # # příklad farmakolog # Kan = read.csv2("data/Kan.csv") summary(Kan) dim(Kan) attach(Kan) plot(y~x);abline(h=0,v=0,lty=3) # funkce nls předpokládá klasikcý zápis regresní funkce # nikoliv zápis "symbolický": bcKan <- nls(y~(x+(625-x)*(1-exp(-b*x/(625-x))))/c,start=list(b=5,c=10)) summary(bcKan) xx <- seq(min(x),max(x),length=101) yy <- predict(bcKan,newdata=list(x=xx)) lines(yy~xx) # podmodel - přímka počátkem summary(cKan <- nls(y~x/c,start=list(c=1))) # zkusme lineární model: summary(lm(y~x-1)) # je náhoda, že t-statistika vyšla identicky? yyp <- predict(cKan,newdata=list(x=xx)) lines(yyp~xx,lty=2) legend("bottomright",legend=c("model","podmodel"),lty=c(1,2)) # test "poměrem věrohodnosti" anova(cKan,bcKan) summary(bcKan) # vyšlo stejně? Porovnej t^2 a F c(summary(bcKan)$coef[1,3]^2,anova(cKan,bcKan)[2,5]) # # byla původní farmakologova parametrizace lepší? bdKan <- nls(y~(x+(625-x)*(1-exp(-b*x/(625-x))))*d,start=list(b=5,d=.1)) summary(bdKan) # porovnání bcKan a bdKan cbind(bcKan=coef(bcKan),bdKan=coef(bdKan)) c(coef(bcKan)[2],"1/"=1/coef(bdKan)[2]) # # zavedení knihovny library(ellipse) # # klasicky výpočtem fce S(theta) pro síť bodů theta fceBC <- function(x,b,c) (x+(625-x)*(1-exp(-b*x/(625-x))))/c fceBD <- function(x,b,d) (x+(625-x)*(1-exp(-b*x/(625-x))))*d SfceBC <- function(b,c) crossprod(y-fceBC(x,b,c)) SfceBD <- function(b,d) crossprod(y-fceBD(x,b,d)) # N <- 101; minD <- 0.15; maxD <- 0.43 #.45 bd <- bc <- matrix(NA,N,N) bb <- seq(0,5,length=N) cc <- seq(1/maxD,1/minD,length=N) dd <- seq(minD,maxD,length=N) for (i in 1:N) for (j in 1:N) bc[i,j] <- SfceBC(bb[i],cc[j]) for (i in 1:N) for (j in 1:N) bd[i,j] <- SfceBD(bb[i],dd[j]) # opar <- par(mfrow = c(1,2), oma = c(1.1, 0, 1.1, 0), las = 1,cex.lab=1.5) contour(bb,cc,bc,levels=deviance(bcKan)*(1+2/30*qf(.95,2,30)), drawlabels=FALSE,lty=1,xlab=expression(beta), ylab=expression(gamma),col="red") points(coef(bcKan)[1],coef(bcKan)[2],pch="+") lines(ellipse(bcKan),lty=2,col="blue") contour(bb,dd,bd,levels=deviance(bdKan)*(1+2/30*qf(.95,2,30)), drawlabels=FALSE,lty=1,xlab=expression(beta), ylab=expression(delta),col="red") points(coef(bdKan)[1],coef(bdKan)[2],pch="+") lines(ellipse(bdKan),lty=2,col="blue") par(opar) # konfidenční intervaly pomocí profilů confint(bcKan) confint(bdKan) 1/confint(bdKan)[2,] # konf. intervaly pomocí Waldova testu (pm = summary(bcKan)$coefficients[,2]*qt(0.975,summary(bcKan)$df[2])) cbind(coef(bcKan)-pm,coef(bcKan)+pm) # profilové diagramy # plot(profile(bcKan,1,alphamax=0.02)) # porovnání s waldovskými intervaly alpha = c(0.01,0.05,0.10,0.20,0.50) (tStar = qt(1-alpha/2,summary(bcKan)$df[2])) (pm = tStar*summary(bcKan)$coeff[1,2]) sort(c(coef(bcKan)[1]-pm,coef(bcKan)[1],coef(bcKan)[1]+pm)) # plot(profile(bcKan,2,alphamax=0.02)) # profilový diagram pro nový parametr plot(profile(bdKan,2,alphamax=0.02)) # současně trojice profilů opar <- par(mfrow = c(1,3), oma = c(1.1, 0, 1.1, 0), las = 1,cex.lab=1.5) plot(profile(bcKan,1,alphamax=0.02)) plot(profile(bcKan,2,alphamax=0.02)) plot(profile(bdKan,2,alphamax=0.02)) par(opar) #