R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Fit a (constrained) piecewise linear regression spline

plrs	R Documentation

Fit a (constrained) piecewise linear regression spline

Description

The function fits a piecewise linear regression spline to explain gene expression by the segmented DNA copy number. The called copy number values are used as a template for model building.

Usage

plrs(expr, cghseg, cghcall=NULL, probloss = NULL, probnorm = NULL,
probgain = NULL, probamp = NULL, knots = NULL, continuous = FALSE,
constr = TRUE, constr.slopes = 2, constr.intercepts = TRUE,
min.obs = 3, discard.obs = TRUE)

Arguments

`expr`	Vector of gene expression values
`cghseg`	Vector of segmented copy number values
`cghcall`	Vector of called copy number values. If not provided, we are reduced to a simple linear model.
`probloss`	Vector of call probabilities associated with state "loss". Default is `NULL`.
`probnorm`	Vector of call probabilities associated with state "normal". Default is `NULL`.
`probgain`	Vector of call probabilities associated with state "gain". Default is `NULL`.
`probamp`	Vector of call probabilities associated with state "amplification". Default is `NULL`.
`knots`	knots or change points. If `NULL` (default), there are estimated. See details.
`continuous`	Logical, whether the model is continuous (no jump) or not.
`constr`	Logical, whether the model is constrained or not. (this has been implemented to turn on and off easily the constraints)
`constr.slopes`	Type of non-negativity constraints applied on slopes. Either 1 or 2 (default). See details.
`constr.intercepts`	If `TRUE` (default) jumps from state to state are also constrained to be non-negative
`min.obs`	See `modify.conf`
`discard.obs`	See `modify.conf`

Details

If cghcall=NULL, discrete copy number values are omitted, which results in fitting a simple linear model.

If constr.slopes=1, all slopes are constrained to be non-negative. If constr.slopes=2, the slope associated with state "normal" is constrained to be non-negative and all others are forced to be at least equal to the latter.

Two methods are implemented for the estimation of knots. If call probabilities are provided, a knot is determined so that the sum of (the two adjacent) states membership probabilities is maximized. Otherwise, this is defined as the midpoint of the interval between the two consecutive states.

The constrained least squares problem is solved using function solve.QP of package quadprog.

Value

An object of class plrs-class

Author(s)

Gwenael G.R. Leday g.g.r.leday@vu.nl

Examples


# Simulate data
sim <- plrs.sim(n=80, states=4, sigma=0.5)

# Fit a model 
model <- plrs(expr=sim$expr, cghseg=sim$seg, cghcall=sim$cal)
model

# Methods
coef(model)
effects(model)
fitted(model)
knots(model)
model.matrix(model)
plot(model)
predict(model, newcghseg=seq(0,5, length.out=100))   
residuals(model)
summary(model)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(plrs)
Loading required package: Biobase
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
    get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
    rbind, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/plrs/plrs.Rd_%03d_medium.png", width=480, height=480)
> ### Name: plrs
> ### Title: Fit a (constrained) piecewise linear regression spline
> ### Aliases: plrs plrs,ANY
> 
> ### ** Examples
> 
> 
> # Simulate data
> sim <- plrs.sim(n=80, states=4, sigma=0.5)
> 
> # Fit a model 
> model <- plrs(expr=sim$expr, cghseg=sim$seg, cghcall=sim$cal)
> model

Object of class "plrs"

Spline coefficients:
theta0.0 theta0.1 theta1.0 theta1.1 theta2.0 theta2.1 theta3.0 theta3.1 
 0.85377  1.17563  0.00000 -1.17563  0.00000  2.91956  0.00000 -1.71605 

Model is constrained:
constr.slopes = 2 
constr.intercepts = TRUE 

> 
> # Methods
> coef(model)
     theta0.0      theta0.1      theta1.0      theta1.1      theta2.0 
 8.537695e-01  1.175631e+00 -2.126995e-16 -1.175631e+00  0.000000e+00 
     theta2.1      theta3.0      theta3.1 
 2.919560e+00 -5.760185e-18 -1.716047e+00 
> effects(model)
       n.obs I.effect S.effect
Loss      20        0    1.176
Normal    20       NA       NA
Gain      20        0    2.920
Ampl.     20        0   -1.716
> fitted(model)
 [1] 1.2737777 0.9460161 1.0900650 1.0877888 1.2076775 1.4071111 1.3847265
 [8] 1.3548981 1.2100264 0.9920386 1.3898524 1.3554946 1.3593043 1.4273300
[15] 1.4383340 1.0779741 1.0837080 1.0853677 1.2268266 1.1254473 1.4414155
[22] 1.4414155 1.4414155 1.4414155 1.4414155 1.4414155 1.4414155 1.4414155
[29] 1.4414155 1.4414155 1.4414155 1.4414155 1.4414155 1.4414155 1.4414155
[36] 1.4414155 1.4414155 1.4414155 1.4414155 1.4414155 1.8732742 2.3236651
[43] 2.1793706 1.4861634 2.6544722 2.4921404 1.5595130 2.0121203 1.6478257
[50] 2.0536496 1.7258645 2.0026541 2.5877203 2.1304619 1.8010788 2.7318525
[57] 2.6423003 1.6414018 2.8419705 2.5431616 3.2579310 2.9829918 3.1373786
[64] 3.1773510 3.5164857 3.2073345 3.0223833 3.0576473 3.2917300 2.9855856
[71] 3.1709777 3.3127617 3.2472772 3.3117200 3.2667284 3.2787483 3.2082847
[78] 3.4834867 3.1781228 3.4469795
> knots(model)
[1] 0.4998559 0.9942567 1.5081741
> model.matrix(model)
      theta0.0   theta0.1 theta1.0    theta1.1 theta2.0   theta2.1 theta3.0
 [1,]        1 0.35726194        0 0.000000000        0 0.00000000        0
 [2,]        1 0.07846556        0 0.000000000        0 0.00000000        0
 [3,]        1 0.20099460        0 0.000000000        0 0.00000000        0
 [4,]        1 0.19905844        0 0.000000000        0 0.00000000        0
 [5,]        1 0.30103661        0 0.000000000        0 0.00000000        0
 [6,]        1 0.47067624        0 0.000000000        0 0.00000000        0
 [7,]        1 0.45163580        0 0.000000000        0 0.00000000        0
 [8,]        1 0.42626356        0 0.000000000        0 0.00000000        0
 [9,]        1 0.30303464        0 0.000000000        0 0.00000000        0
[10,]        1 0.11761268        0 0.000000000        0 0.00000000        0
[11,]        1 0.45599591        0 0.000000000        0 0.00000000        0
[12,]        1 0.42677089        0 0.000000000        0 0.00000000        0
[13,]        1 0.43001144        0 0.000000000        0 0.00000000        0
[14,]        1 0.48787464        0 0.000000000        0 0.00000000        0
[15,]        1 0.49723473        0 0.000000000        0 0.00000000        0
[16,]        1 0.19070996        0 0.000000000        0 0.00000000        0
[17,]        1 0.19558731        0 0.000000000        0 0.00000000        0
[18,]        1 0.19699902        0 0.000000000        0 0.00000000        0
[19,]        1 0.31732498        0 0.000000000        0 0.00000000        0
[20,]        1 0.23109101        0 0.000000000        0 0.00000000        0
[21,]        1 0.74750367        1 0.247647792        0 0.00000000        0
[22,]        1 0.87425018        1 0.374394305        0 0.00000000        0
[23,]        1 0.79720410        1 0.297348229        0 0.00000000        0
[24,]        1 0.96921378        1 0.469357900        0 0.00000000        0
[25,]        1 0.50247702        1 0.002621147        0 0.00000000        0
[26,]        1 0.66434349        1 0.164487618        0 0.00000000        0
[27,]        1 0.60211738        1 0.102261500        0 0.00000000        0
[28,]        1 0.72422358        1 0.224367705        0 0.00000000        0
[29,]        1 0.97892972        1 0.479073841        0 0.00000000        0
[30,]        1 0.51386066        1 0.014004782        0 0.00000000        0
[31,]        1 0.96812578        1 0.468269909        0 0.00000000        0
[32,]        1 0.59523346        1 0.095377589        0 0.00000000        0
[33,]        1 0.62636926        1 0.126513385        0 0.00000000        0
[34,]        1 0.61920093        1 0.119345050        0 0.00000000        0
[35,]        1 0.67509202        1 0.175236145        0 0.00000000        0
[36,]        1 0.64423744        1 0.144381562        0 0.00000000        0
[37,]        1 0.80503607        1 0.305180195        0 0.00000000        0
[38,]        1 0.63699710        1 0.137141223        0 0.00000000        0
[39,]        1 0.93877712        1 0.438921240        0 0.00000000        0
[40,]        1 0.90155967        1 0.401703798        0 0.00000000        0
[41,]        1 1.14217576        1 0.642319887        1 0.14791911        0
[42,]        1 1.29644247        1 0.796586598        1 0.30218582        0
[43,]        1 1.24701911        1 0.747163237        1 0.25276246        0
[44,]        1 1.00958359        1 0.509727710        1 0.01532693        0
[45,]        1 1.40974966        1 0.909893783        1 0.41549301        0
[46,]        1 1.35414821        1 0.854292337        1 0.35989156        0
[47,]        1 1.03470708        1 0.534851203        1 0.04045043        0
[48,]        1 1.18973297        1 0.689877090        1 0.19547631        0
[49,]        1 1.06495572        1 0.565099841        1 0.07069907        0
[50,]        1 1.20395747        1 0.704101590        1 0.20970081        0
[51,]        1 1.09168537        1 0.591829498        1 0.09742872        0
[52,]        1 1.18649063        1 0.686634758        1 0.19223398        0
[53,]        1 1.38688598        1 0.887030102        1 0.39262933        0
[54,]        1 1.23026702        1 0.730411145        1 0.23601037        0
[55,]        1 1.11744758        1 0.617591703        1 0.12319093        0
[56,]        1 1.43625377        1 0.936397898        1 0.44199712        0
[57,]        1 1.40558058        1 0.905724700        1 0.41132392        0
[58,]        1 1.06275543        1 0.562899549        1 0.06849877        0
[59,]        1 1.47397110        1 0.974115226        1 0.47971445        0
[60,]        1 1.37162384        1 0.871767966        1 0.37736719        0
[61,]        1 1.77082426        1 1.270968382        1 0.77656761        1
[62,]        1 1.54237701        1 1.042521139        1 0.54812036        1
[63,]        1 1.67065721        1 1.170801336        1 0.67640056        1
[64,]        1 1.70387032        1 1.204014446        1 0.70961367        1
[65,]        1 1.98565763        1 1.485801753        1 0.99140098        1
[66,]        1 1.72878361        1 1.228927738        1 0.73452696        1
[67,]        1 1.57510751        1 1.075251639        1 0.58085086        1
[68,]        1 1.60440838        1 1.104552508        1 0.61015173        1
[69,]        1 1.79890790        1 1.299052021        1 0.80465125        1
[70,]        1 1.54453225        1 1.044676377        1 0.55027560        1
[71,]        1 1.69857474        1 1.198718862        1 0.70431809        1
[72,]        1 1.81638317        1 1.316527291        1 0.82212652        1
[73,]        1 1.76197202        1 1.262116147        1 0.76771537        1
[74,]        1 1.81551764        1 1.315661763        1 0.82126099        1
[75,]        1 1.77813407        1 1.278278199        1 0.78387742        1
[76,]        1 1.78812144        1 1.288265560        1 0.79386478        1
[77,]        1 1.72957316        1 1.229717284        1 0.73531651        1
[78,]        1 1.95823875        1 1.458382879        1 0.96398210        1
[79,]        1 1.70451156        1 1.204655682        1 0.71025491        1
[80,]        1 1.92790487        1 1.428048996        1 0.93364822        1
        theta3.1
 [1,] 0.00000000
 [2,] 0.00000000
 [3,] 0.00000000
 [4,] 0.00000000
 [5,] 0.00000000
 [6,] 0.00000000
 [7,] 0.00000000
 [8,] 0.00000000
 [9,] 0.00000000
[10,] 0.00000000
[11,] 0.00000000
[12,] 0.00000000
[13,] 0.00000000
[14,] 0.00000000
[15,] 0.00000000
[16,] 0.00000000
[17,] 0.00000000
[18,] 0.00000000
[19,] 0.00000000
[20,] 0.00000000
[21,] 0.00000000
[22,] 0.00000000
[23,] 0.00000000
[24,] 0.00000000
[25,] 0.00000000
[26,] 0.00000000
[27,] 0.00000000
[28,] 0.00000000
[29,] 0.00000000
[30,] 0.00000000
[31,] 0.00000000
[32,] 0.00000000
[33,] 0.00000000
[34,] 0.00000000
[35,] 0.00000000
[36,] 0.00000000
[37,] 0.00000000
[38,] 0.00000000
[39,] 0.00000000
[40,] 0.00000000
[41,] 0.00000000
[42,] 0.00000000
[43,] 0.00000000
[44,] 0.00000000
[45,] 0.00000000
[46,] 0.00000000
[47,] 0.00000000
[48,] 0.00000000
[49,] 0.00000000
[50,] 0.00000000
[51,] 0.00000000
[52,] 0.00000000
[53,] 0.00000000
[54,] 0.00000000
[55,] 0.00000000
[56,] 0.00000000
[57,] 0.00000000
[58,] 0.00000000
[59,] 0.00000000
[60,] 0.00000000
[61,] 0.26265020
[62,] 0.03420296
[63,] 0.16248315
[64,] 0.19569626
[65,] 0.47748357
[66,] 0.22060956
[67,] 0.06693346
[68,] 0.09623433
[69,] 0.29073384
[70,] 0.03635819
[71,] 0.19040068
[72,] 0.30820911
[73,] 0.25379796
[74,] 0.30734358
[75,] 0.26996002
[76,] 0.27994738
[77,] 0.22139910
[78,] 0.45006470
[79,] 0.19633750
[80,] 0.41973081
> plot(model)
> predict(model, newcghseg=seq(0,5, length.out=100))   
  [1] 0.8537695 0.9131448 0.9725201 1.0318954 1.0912707 1.1506460 1.2100213
  [8] 1.2693966 1.3287719 1.3881472 1.4414155 1.4414155 1.4414155 1.4414155
 [15] 1.4414155 1.4414155 1.4414155 1.4414155 1.4414155 1.4414155 1.4876741
 [22] 1.6351266 1.7825791 1.9300316 2.0774841 2.2249366 2.3723891 2.5198416
 [29] 2.6672942 2.8147467 2.9502255 3.0110090 3.0717925 3.1325760 3.1933594
 [36] 3.2541429 3.3149264 3.3757099 3.4364934 3.4972768 3.5580603 3.6188438
 [43] 3.6796273 3.7404108 3.8011942 3.8619777 3.9227612 3.9835447 4.0443281
 [50] 4.1051116 4.1658951 4.2266786 4.2874621 4.3482455 4.4090290 4.4698125
 [57] 4.5305960 4.5913795 4.6521629 4.7129464 4.7737299 4.8345134 4.8952968
 [64] 4.9560803 5.0168638 5.0776473 5.1384308 5.1992142 5.2599977 5.3207812
 [71] 5.3815647 5.4423482 5.5031316 5.5639151 5.6246986 5.6854821 5.7462655
 [78] 5.8070490 5.8678325 5.9286160 5.9893995 6.0501829 6.1109664 6.1717499
 [85] 6.2325334 6.2933169 6.3541003 6.4148838 6.4756673 6.5364508 6.5972342
 [92] 6.6580177 6.7188012 6.7795847 6.8403682 6.9011516 6.9619351 7.0227186
 [99] 7.0835021 7.1442856
> residuals(model)
 [1]  0.73673381 -0.98959433 -0.01979445  0.71549996 -0.12659880 -0.39408500
 [7] -0.79422586  0.31675618  0.40860538  0.12688818  0.66178626  0.65809893
[13]  0.26359331 -0.44283671 -0.22855407 -0.51870572 -0.35707817  0.40711297
[19]  0.34955033  0.06568142  0.45316310  0.17388659 -0.53007494  0.45147842
[25] -0.29681926 -0.04894283  1.02252505 -0.77859493 -1.00817760 -0.29205339
[31]  0.02640543  0.85481692  0.29444597  0.27775816 -0.03783276 -0.46155957
[37] -1.43096769  1.11165124 -0.64280115  0.22488252  0.12000587  0.07550985
[43] -0.02479635  0.12778712  0.10444892  0.31478852 -0.19857829  0.07045400
[49] -0.34711759 -0.30475855  0.31456975  0.25450887 -0.32512775  0.66996878
[55] -0.08630388  0.27894659  0.68922694 -0.67415445  0.04285270 -0.58195997
[61]  0.14556448  0.53813079  0.05830690 -0.19938995 -0.23086908 -0.99420998
[67] -0.15266880 -0.59320279  0.90669070 -0.33507574 -0.23000294  0.26682042
[73] -0.23099177 -0.66209270  0.38371807  0.11520423  0.40927043 -0.03761709
[79] -0.36746736  0.48758820
> summary(model)

Object of class "plrs"

Spline coefficients:
theta0.0 theta0.1 theta1.0 theta1.1 theta2.0 theta2.1 theta3.0 theta3.1 
 0.85377  1.17563  0.00000 -1.17563  0.00000  2.91956  0.00000 -1.71605 

Model is constrained:
constr.slopes = 2 
constr.intercepts = TRUE 

Configuration:
       n.obs I.effect S.effect
Loss      20        0    1.176
Normal    20       NA       NA
Gain      20        0    2.920
Ampl.     20        0   -1.716
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>