R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Laajennuspaketti drc: EC50-arvon määrittäminen

Effective concentration 50% (EC50)-arvolla mitataan esimerkiksi lääkeaineiden vaikutusta. EC50-arvo on lääkeaineen konsentraatio, jolla saadaan aikaan puolet vaikutuksen maksimiarvosta. EC50-arvon lähisukulainen on effective dose (ED50), joka taas vastaa konsentraatiota, jolla saadaan vaste 50 %:ssa altistetuista kohteista, esimerkiksi potilaista. Yleensä ED50-arvoa käytetään hyödyllisistä, tavoitelluista terapeuttisista vaikutuksista. Jos tutkitaan aineiden myrkyllisyyksiä, käytetään arvoa LD50, joka vastaa konsentraatiota, jossa 50 % altistetuista kohteista kuolee.

EC50-arvon määrittämiseen on useita erilaisia tapoja. Olen aiemmin kirjoittanut epäparametrisesta Kärberin menetelmästä. Usein sovellettu probit-malli voidaan sovittaa funktiolla glm(), mutta laajennuspaketti drc tarjoaa laajemman valikoiman erilaisia annosvasteaineistojen analysoimiseen soveltuvia menetelmiä.

Seuraavassa käytetään esimerkkinä aineistoa, jossa soluviljelmä on altistettu myrkylliselle aineelle. Määritetään EC50/ED50-arvo, jolla puolet soluista selviää käsittelystä.

Muodostetaan ensin esimerkkiaineisto objektiin dat:

dat<-structure(list(conc = c(0L, 0L, 0L, 50L, 50L, 50L, 100L, 100L, 
100L, 200L, 200L, 200L, 300L, 300L, 300L, 400L, 400L, 400L, 500L, 
500L, 500L, 750L, 750L, 750L, 1000L, 1000L, 1000L, 5000L, 5000L, 
5000L, 10000L, 10000L, 10000L), resp = c(21250L, 19750L, 20250L, 
18900L, 20850L, 16850L, 16800L, 19350L, 19100L, 17550L, 18650L, 
17550L, 17350L, 19750L, 18200L, 17150L, 17350L, 18700L, 14500L, 
13900L, 15650L, 12600L, 11150L, 11000L, 8025L, 8250L, 9150L, 
159L, 111L, 129L, 12L, 0L, 0L)), .Names = c("conc", "resp"), class = "data.frame", row.names = c(NA, 
-33L))
 
print(dat)
 
    conc  resp
1      0 21250
2      0 19750
3      0 20250
4     50 18900
5     50 20850
6     50 16850
7    100 16800
8    100 19350
9    100 19100
10   200 17550
11   200 18650
12   200 17550
13   300 17350
14   300 19750
15   300 18200
16   400 17150
17   400 17350
18   400 18700
19   500 14500
20   500 13900
21   500 15650
22   750 12600
23   750 11150
24   750 11000
25  1000  8025
26  1000  8250
27  1000  9150
28  5000   159
29  5000   111
30  5000   129
31 10000    12
32 10000     0
33 10000     0

Aineistossa on kaksi saraketta, joista resp on käsittelyn jälkeen laskettu solulukumäärä ja conc on viljelmään lisätyn myrkyllisen aineen konsentraatio.

Aineistoon voidaan sovittaa vaikkapa neljä erilaista logistista mallia, joista lopulliseksi valitaan malli, jonka AIC-arvo on kaikkein pienin. Paketti drc tarjoaa monia muitakin mahdollisia malleja kuin logistisen mallin. Mukana ovat mm. log-logistinen- ja Weibull-mallit. Käytettävä malli määritellään sen lyhenteen perusteella. Esimerkiksi logistinen malli, riippuen sen parametrien määrästä, on nimetty joko L2, L3, L4 tai L5. Seuraavassa sovitetaan mainitut logistiset mallit aineistoon, ja lasketaan niiden AIC-arvot sekä tallennetaan paras malli objektiin best. Lopuksi parasta mallia käytetään EC50-arvon laskemiseen:

library(drc)
 
fit2<-fit<-drm(resp~conc, fct=LL.2(), data=dat)
fit3<-fit<-drm(resp~conc, fct=LL.3(), data=dat)
fit4<-fit<-drm(resp~conc, fct=LL.4(), data=dat)
fit5<-fit<-drm(resp~conc, fct=LL.5(), data=dat)
 
AIC(fit2, fit3, fit4, fit5)
 
#     df      AIC
#fit2  3 734.8355
#fit3  4 559.8222  # Paras malli
#fit4  5 561.7956
#fit5  6 563.7944
 
best<-switch(which(abs(diff(AIC(fit2, fit3, fit4, fit5)$AIC))<2)[1]+1, best=fit1, best=fit2, best=fit3, best=fit4)       
 
ed<-ED(best, c(50))
 
#Estimated effective doses
#
#       Estimate Std. Error
#e:1:50   895.46      37.06

Tässä aineistossa EC50-arvo on noin 895. Kun myrkylliseen aineen konsentraatio on noin 895, se tappaa suunnilleen puolet soluviljelmän soluista.

Eri mallien antamat tulokset voidaan myös visualisoida esimerkiksi seuraavasti:

par(mfrow=c(2,2), mar=c(4,5,3,1))
plot(fit2, main="2-parameter logistic model")
plot(fit3, main="3-parameter logistic model")
plot(fit4, main="4-parameter logistic model")
plot(fit5, main="5-parameter logistic model")

Muodostuva kuva on seuraavanlainen:

drc1

Kuvasta on helppo nähdä, että kunhan logistisessa mallissa on vähintään kolme parametria, ovat sovitteet keskenään varsin samannäköisiä.

Parhaan mallin antama tulos voidaan visualisoida myös virhetermeineen tarkemmin esimerkiksi seuraavalla tavalla. Samaan kuvaan voidaan merkitä myös EC50-arvo ja sen keskihajonta:

plot(best, lwd=2, pch=16, col="black", main=paste("Model: ", best$fct$name, ", AIC: ", round(AIC(best), 2), sep=""), xlab="", ylab="", cex=0.75)
for(i in 1:length(unique(dat$conc))) {
   if(i==1) {
      lines(x=c(unique(dat$conc)[i]+10, unique(dat$conc)[i]+10), y=c(aggregate(dat$resp, list(dat$conc), mean)$x[i]-aggregate(dat$resp, list(dat$conc), sd)$x[i]*1.96, aggregate(dat$resp, list(dat$conc), mean)$x[i]+aggregate(dat$resp, list(dat$conc), sd)$x[i]*1.96))
   } else {
      lines(x=c(unique(dat$conc)[i], unique(dat$conc)[i]), y=c(aggregate(dat$resp, list(dat$conc), mean)$x[i]-aggregate(dat$resp, list(dat$conc), sd)$x[i]*1.96, aggregate(dat$resp, list(dat$conc), mean)$x[i]+aggregate(dat$resp, list(dat$conc), sd)$x[i]*1.96))
   }
}
lines(x=c(ed[1], ed[1]), y=c(-diff(range(dat$conc))*1.04, predict(best, data.frame(conc=ed[1]))), col="grey75") 
lines(x=c(ed[1]+ed[2], ed[1]+ed[2]), y=c(-diff(range(dat$conc))*1.04, predict(best, data.frame(conc=ed[1]+ed[2]))), col="grey75", lty=2) 
lines(x=c(ed[1]-ed[2], ed[1]-ed[2]), y=c(-diff(range(dat$conc))*1.04, predict(best, data.frame(conc=ed[1]-ed[2]))), col="grey75", lty=2) 
plot(best, lwd=2, pch=16, col="red", type="none", add=TRUE)
mtext("Concentration", side=1, las=1, font=2, line=3)
mtext("Response", side=2, at=max(dat$resp)*1.1, las=1, font=2)
legend(x="topright", legend=paste("EC50 = ", round(ed[1], 1), sep=""), bty="n", cex=1) 
box()

Esimerkkidatasta piirretty kuva näyttää seuraavalta:

drc2

Paketti on drc on varsin monikäyttöinen, sillä sillä voidaan sovittaa myös Michaelis-Mentenin malleja mittausdataan. Lisäksi paketin funktioilla voidaan sovittaa ns. NEC (no effect concentration) -malleja, joissa oletetaan olevan jokin konsentraatioraja, jonka alle jäävillä arvoilla ei vastekaan eroa peruslinjasta tai kontrolleista millään tavalla. Pakettiin kannattaa ehdottomasti tutustua, jos funktion glm() tarjoamat mahdollisuudet eivät riitä täyttämään annovasteanalyysien tarpeita!

Tags: ,