Desktop version

Home arrow Economics

  • Increase font
  • Decrease font


<<   CONTENTS   >>

Analysis of One Gene Using the gls Function

The main function in the package, fitJM(), fits a joint model for the gene expression and the bioactivity data. In this section we discuss the implementation in R for an analysis of a particular gene using the gls() function. For illustration we use the data for the gene FGFBP1 in the EGFR dataset. The object data contains the variables Responses (gene expression and bioassay data), myCov (the fingerprint feature) sampleID and respIndex (endpoint identification: 1—gene expression, 2—bioassay data). A partial print of the data is given below.

>head(data)

Responses myCov samplelD respIndex 1 11.08821 011

  • 2 10.18648 121
  • 3 11.22024 031
  • 4 10.74484 141
  • 5 11.08644 151
  • 6 10.87207 1 6 1

The joint model specified in (16.1) can be fitted using the gls() function in the following way:

>f1<- gls(Responses ~ myCov*as.factor(respIndex),

correlation = corSymm(form = ~ respIndex| sampleID), weights = varIdent(form=~ 1|respIndex), data = data, method = "ML")

Parameter estimates for the chemical structure effects, a = 0.7872531 and в = -1.727113 + 0.787253 = -0.93986, are shown in the panel below. These values are reported in the list of top K genes and in Table 16.4.

summary(f1)

Coefficients:

Value Std.Error t-value p-value (Intercept) 10.116852 0.1776370 56.95240 0e+00

myCov 0.787253 0.2240557 3.51365 8e-04

as.factor(respIndex)2 -3.926083 0.3043475 -12.90000 0e+00

myCov:as.factor(respIndex)2 -1.727113 0.3838773 -4.49913 0e+00

The adjusted association, p = -0.781, can be found in the correlation structure panel:

Correlation Structure: General Formula: ~respIndex | sampleID Parameter estimate(s):

Correlation:

  • 1
  • 2 -0.781

To test the null hypothesis H0 : p = 0, we fit a reduced model in which the covariance between gene expression and the bioactivity data is zero in the following way:

>f2 <- gls(Responses ~ myCov*as.factor(respIndex), weights = varIdent(form=~ 1|respIndex), data = data, method="ML")

The likelihood ratio test is used to test the null hypothesis. For the gene FGFBP1, the null hypothesis is rejected:

>anova(f1,f2) Model df AIC BIC logLik Test L.Ratio p-value f1 17 98.77017 114.5096 -42.38508

f2 26 129.77289 143.2639 -58.88644 1 vs 2 33.00272 <.0001

 
<<   CONTENTS   >>

Related topics