python - RPy 中 R 的 pdIndent 函数

Question

我正在努力将 lmeSplines tutorial 的代码翻译为 RPy。

我现在陷入了以下行:

fit1s <- lme(y ~ time, data=smSplineEx1,random=list(all=pdIdent(~Zt - 1)))

我之前曾使用过nlme.lme，并且以下工作正常:

from rpy2.robjects.packages import importr
nlme = importr('nlme')    
nlme.lme(r.formula('y ~ time'), data=some_data, random=r.formula('~1|ID'))

但是这还有另一个随机分配。我想知道如何翻译这一点并将其放入我的 RPy 代码中 list(all=pdIdent(~Zt - 1))。

(预处理的)示例数据 smSplineEx1 的结构如下所示(Zt.* 高达 98):

    time         y   y.true all          Zt.1          Zt.2          Zt.3
1      1  5.797149 4.235263   1  1.168560e+00  2.071261e+00  2.944953e+00
2      2  5.469222 4.461302   1  1.487859e-01  1.072013e+00  1.948857e+00
3      3  4.567237 4.678477   1 -5.449190e-02  7.276623e-02  9.527613e-01
4      4  3.645763 4.887137   1 -5.364552e-02 -1.359115e-01 -4.333438e-02
5      5  5.094126 5.087615   1 -5.279913e-02 -1.337708e-01 -2.506194e-01
6      6  4.636121 5.280233   1 -5.195275e-02 -1.316300e-01 -2.466158e-01
7      7  5.501538 5.465298   1 -5.110637e-02 -1.294892e-01 -2.426123e-01
8      8  5.011509 5.643106   1 -5.025998e-02 -1.273485e-01 -2.386087e-01
9      9  6.114037 5.813942   1 -4.941360e-02 -1.252077e-01 -2.346052e-01
10    10  5.696472 5.978080   1 -4.856722e-02 -1.230670e-01 -2.306016e-01
11    11  6.615363 6.135781   1 -4.772083e-02 -1.209262e-01 -2.265980e-01
12    12  8.002526 6.287300   1 -4.687445e-02 -1.187854e-01 -2.225945e-01
13    13  6.887444 6.432877   1 -4.602807e-02 -1.166447e-01 -2.185909e-01
14    14  6.319205 6.572746   1 -4.518168e-02 -1.145039e-01 -2.145874e-01
15    15  6.482771 6.707130   1 -4.433530e-02 -1.123632e-01 -2.105838e-01
16    16  7.938015 6.836245   1 -4.348892e-02 -1.102224e-01 -2.065802e-01
17    17  7.585533 6.960298   1 -4.264253e-02 -1.080816e-01 -2.025767e-01
18    18  7.560287 7.079486   1 -4.179615e-02 -1.059409e-01 -1.985731e-01
19    19  7.571020 7.194001   1 -4.094977e-02 -1.038001e-01 -1.945696e-01
20    20  8.922418 7.304026   1 -4.010338e-02 -1.016594e-01 -1.905660e-01
21    21  8.241394 7.409737   1 -3.925700e-02 -9.951861e-02 -1.865625e-01
22    22  7.447076 7.511303   1 -3.841062e-02 -9.737785e-02 -1.825589e-01
23    23  7.317292 7.608886   1 -3.756423e-02 -9.523709e-02 -1.785553e-01
24    24  7.077333 7.702643   1 -3.671785e-02 -9.309633e-02 -1.745518e-01
25    25  8.268601 7.792723   1 -3.587147e-02 -9.095557e-02 -1.705482e-01
26    26  8.216013 7.879272   1 -3.502508e-02 -8.881481e-02 -1.665447e-01
27    27  8.968495 7.962427   1 -3.417870e-02 -8.667405e-02 -1.625411e-01
28    28  9.085605 8.042321   1 -3.333232e-02 -8.453329e-02 -1.585375e-01
29    29  9.002575 8.119083   1 -3.248593e-02 -8.239253e-02 -1.545340e-01
30    30  8.763187 8.192835   1 -3.163955e-02 -8.025177e-02 -1.505304e-01
31    31  8.936370 8.263695   1 -3.079317e-02 -7.811101e-02 -1.465269e-01
32    32  9.033403 8.331776   1 -2.994678e-02 -7.597025e-02 -1.425233e-01
33    33  8.248328 8.397188   1 -2.910040e-02 -7.382949e-02 -1.385198e-01
34    34  5.961721 8.460035   1 -2.825402e-02 -7.168873e-02 -1.345162e-01
35    35  8.400489 8.520418   1 -2.740763e-02 -6.954797e-02 -1.305126e-01
36    36  6.855125 8.578433   1 -2.656125e-02 -6.740721e-02 -1.265091e-01
37    37  9.798931 8.634174   1 -2.571487e-02 -6.526645e-02 -1.225055e-01
38    38  8.862758 8.687729   1 -2.486848e-02 -6.312569e-02 -1.185020e-01
39    39  7.282970 8.739184   1 -2.402210e-02 -6.098493e-02 -1.144984e-01
40    40  7.484208 8.788621   1 -2.317572e-02 -5.884417e-02 -1.104949e-01
41    41  8.404670 8.836120   1 -2.232933e-02 -5.670341e-02 -1.064913e-01
42    42  8.880734 8.881756   1 -2.148295e-02 -5.456265e-02 -1.024877e-01
43    43  8.826189 8.925603   1 -2.063657e-02 -5.242189e-02 -9.848418e-02
44    44  9.827906 8.967731   1 -1.979018e-02 -5.028113e-02 -9.448062e-02
45    45  8.528795 9.008207   1 -1.894380e-02 -4.814037e-02 -9.047706e-02
46    46  9.484073 9.047095   1 -1.809742e-02 -4.599961e-02 -8.647351e-02
47    47  8.911947 9.084459   1 -1.725103e-02 -4.385885e-02 -8.246995e-02
48    48 10.201343 9.120358   1 -1.640465e-02 -4.171809e-02 -7.846639e-02
49    49  8.908016 9.154849   1 -1.555827e-02 -3.957733e-02 -7.446283e-02
50    50  8.202368 9.187988   1 -1.471188e-02 -3.743657e-02 -7.045927e-02
51    51  7.432851 9.219828   1 -1.386550e-02 -3.529581e-02 -6.645572e-02
52    52  8.063268 9.250419   1 -1.301912e-02 -3.315505e-02 -6.245216e-02
53    53 10.155756 9.279810   1 -1.217273e-02 -3.101429e-02 -5.844860e-02
54    54  7.905281 9.308049   1 -1.132635e-02 -2.887353e-02 -5.444504e-02
55    55  9.688337 9.335181   1 -1.047997e-02 -2.673277e-02 -5.044148e-02
56    56  9.437176 9.361249   1 -9.633582e-03 -2.459201e-02 -4.643793e-02
57    57  9.165873 9.386295   1 -8.787198e-03 -2.245125e-02 -4.243437e-02
58    58  9.120195 9.410358   1 -7.940815e-03 -2.031049e-02 -3.843081e-02
59    59  9.955840 9.433479   1 -7.094432e-03 -1.816973e-02 -3.442725e-02
60    60  9.314230 9.455692   1 -6.248048e-03 -1.602897e-02 -3.042369e-02
61    61  9.706852 9.477035   1 -5.401665e-03 -1.388821e-02 -2.642014e-02
62    62  9.615765 9.497541   1 -4.555282e-03 -1.174746e-02 -2.241658e-02
63    63  7.918843 9.517242   1 -3.708898e-03 -9.606695e-03 -1.841302e-02
64    64  9.352892 9.536172   1 -2.862515e-03 -7.465935e-03 -1.440946e-02
65    65  9.722685 9.554359   1 -2.016132e-03 -5.325176e-03 -1.040590e-02
66    66  9.186888 9.571832   1 -1.169748e-03 -3.184416e-03 -6.402346e-03
67    67  8.652299 9.588621   1 -3.233650e-04 -1.043656e-03 -2.398788e-03
68    68  8.681421 9.604751   1  5.230184e-04  1.097104e-03  1.604770e-03
69    69 10.279181 9.620249   1  1.369402e-03  3.237864e-03  5.608328e-03
70    70  9.314963 9.635140   1  2.215785e-03  5.378623e-03  9.611886e-03
71    71  6.897151 9.649446   1  3.062168e-03  7.519383e-03  1.361544e-02
72    72  9.343135 9.663191   1  3.908552e-03  9.660143e-03  1.761900e-02
73    73  9.273135 9.676398   1  4.754935e-03  1.180090e-02  2.162256e-02
74    74 10.041796 9.689086   1  5.601318e-03  1.394166e-02  2.562612e-02
75    75  9.724713 9.701278   1  6.447702e-03  1.608242e-02  2.962968e-02
76    76  8.593517 9.712991   1  7.294085e-03  1.822318e-02  3.363323e-02
77    77  7.401988 9.724244   1  8.140468e-03  2.036394e-02  3.763679e-02
78    78 10.258688 9.735057   1  8.986852e-03  2.250470e-02  4.164035e-02
79    79 10.037192 9.745446   1  9.833235e-03  2.464546e-02  4.564391e-02
80    80  9.637510 9.755427   1  1.067962e-02  2.678622e-02  4.964747e-02
81    81  8.887625 9.765017   1  1.152600e-02  2.892698e-02  5.365102e-02
82    82  9.922013 9.774230   1  1.237239e-02  3.106774e-02  5.765458e-02
83    83 10.466709 9.783083   1  1.321877e-02  3.320850e-02  6.165814e-02
84    84 11.132830 9.791588   1  1.406515e-02  3.534926e-02  6.566170e-02
85    85 10.154038 9.799760   1  1.491154e-02  3.749002e-02  6.966526e-02
86    86 10.433068 9.807612   1  1.575792e-02  3.963078e-02  7.366881e-02
87    87  9.666781 9.815156   1  1.660430e-02  4.177154e-02  7.767237e-02
88    88  9.478004 9.822403   1  1.745069e-02  4.391230e-02  8.167593e-02
89    89 10.002749 9.829367   1  1.829707e-02  4.605306e-02  8.567949e-02
90    90  7.593259 9.836058   1  1.914345e-02  4.819382e-02  8.968305e-02
91    91 10.915754 9.842486   1  1.998984e-02  5.033458e-02  9.368660e-02
92    92  8.855580 9.848662   1  2.083622e-02  5.247534e-02  9.769016e-02
93    93  8.884683 9.854596   1  2.168260e-02  5.461610e-02  1.016937e-01
94    94  9.757451 9.860298   1  2.252899e-02  5.675686e-02  1.056973e-01
95    95 10.222361 9.865775   1  2.337537e-02  5.889762e-02  1.097008e-01
96    96  9.090410 9.871038   1  2.422175e-02  6.103838e-02  1.137044e-01
97    97  8.837872 9.876095   1  2.506814e-02  6.317914e-02  1.177080e-01
98    98  9.413135 9.880953   1  2.591452e-02  6.531990e-02  1.217115e-01
99    99  9.295531 9.885621   1  2.676090e-02  6.746066e-02  1.257151e-01
100  100  9.698118 9.890106   1  2.760729e-02  6.960142e-02  1.297186e-01

Answer 1

您可以输入list(all=pdIdent(~Zt - 1))在 R的全局环境使用 reval()方法:

In [55]:

import rpy2.robjects as ro
import pandas.rpy.common as com
mydata = ro.r['data.frame']
read = ro.r['read.csv']
head = ro.r['head']
summary = ro.r['summary']
library = ro.r['library']
In [56]:

formula = '~ time'
library('lmeSplines')
ro.reval('data(smSplineEx1)')
ro.reval('smSplineEx1$all <- rep(1,nrow(smSplineEx1))')
ro.reval('smSplineEx1$Zt <- smspline(~ time, data=smSplineEx1)')
ro.reval('rnd <- list(all=pdIdent(~Zt - 1))')
#result = ro.r.smspline(formula=ro.r(formula), data=ro.r.smSplineEx1) #notice: data=ro.r.smSplineEx1
result = ro.r.lme(ro.r('y~time'), data=ro.r.smSplineEx1, random=ro.r.rnd)
In [57]:

print com.convert_robj(result.rx('coefficients'))
{'coefficients': {'random': {'all':         Zt1       Zt2       Zt3       Zt4       Zt5       Zt6       Zt7  \
1  0.000509  0.001057  0.001352  0.001184  0.000869  0.000283 -0.000424   

        Zt8       Zt9      Zt10    ...         Zt89      Zt90      Zt91  \
1 -0.001367 -0.002325 -0.003405    ...    -0.001506 -0.001347 -0.000864   

       Zt92      Zt93      Zt94      Zt95      Zt96      Zt97      Zt98  
1 -0.000631 -0.000569 -0.000392 -0.000049  0.000127  0.000114  0.000071  

[1 rows x 98 columns]}, 'fixed': (Intercept)    6.498800
time           0.038723
dtype: float64}}

小心，结果有点变形。基本上它是嵌套字典，无法转换为 pandas.DataFrame .

您可以访问y在smsSplineEx通过ro.r.smSplineEx1.rx('y') ，类似于smsSplineEx1$y正如您在 R 中所做的那样.

现在假设您有 result变量 python ，由生成

result = ro.r.lme(ro.r('y~time'), data=ro.r.smSplineEx1, random=ro.r.rnd)

并且您想使用 R 来绘制它，(而不是使用 matplotlib 来绘制它)，您需要将其分配给 R 中的变量工作区:

ro.R().assign('result', result)

现在有一个名为 result 的变量在R工作区，您可以使用 ro.r.result 访问它.

使用 R 绘制它:

In [17]:

ro.reval('plot(smSplineEx1$time,smSplineEx1$y,pch="o",type="n", \
          main="Spline fits: lme(y ~ time, random=list(all=pdIdent(~Zt-1)))", \
          xlab="time",ylab="y")')
Out[17]:
rpy2.rinterface.NULL
In [21]:

ro.reval('lines(smSplineEx1$time, fitted(result),col=2)')
Out[21]:
rpy2.rinterface.NULL

或者您可以在 R 中执行所有操作:

ro.reval('result <- lme(y ~ time, data=smSplineEx1,random=list(all=pdIdent(~Zt - 1)))')
ro.reval('plot(smSplineEx1$time,smSplineEx1$y,pch="o",type="n", \
         main="Spline fits: lme(y ~ time, random=list(all=pdIdent(~Zt-1)))", \
         xlab="time",ylab="y")')
ro.reval('lines(smSplineEx1$time, fitted(result),col=2)')

并访问 R变量使用:ro.r.smSplineEx1.rx2('time')或ro.r.result

编辑

注意一些R对象无法转换为 pandas.dataFrame由于数据结构的混合而保持原样:

In [62]:

ro.r["smSplineEx1"]
Out[62]:
<DataFrame - Python:0x108525518 / R:0x109e5da38>
[FloatVe..., FloatVe..., FloatVe..., FloatVe..., Matrix]
  time: <class 'rpy2.robjects.vectors.FloatVector'>
  <FloatVector - Python:0x10807e518 / R:0x1022599e0>
[1.000000, 2.000000, 3.000000, ..., 98.000000, 99.000000, 100.000000]
  y: <class 'rpy2.robjects.vectors.FloatVector'>
  <FloatVector - Python:0x108525a70 / R:0x102259d30>
[5.797149, 5.469222, 4.567237, ..., 9.413135, 9.295531, 9.698118]
  y.true: <class 'rpy2.robjects.vectors.FloatVector'>
  <FloatVector - Python:0x1085257a0 / R:0x10225dfb0>
[4.235263, 4.461302, 4.678477, ..., 9.880953, 9.885621, 9.890106]
  all: <class 'rpy2.robjects.vectors.FloatVector'>
  <FloatVector - Python:0x1085258c0 / R:0x10225e300>
[1.000000, 1.000000, 1.000000, ..., 1.000000, 1.000000, 1.000000]
  Zt: <class 'rpy2.robjects.vectors.Matrix'>
  <Matrix - Python:0x108525908 / R:0x103e8ba00>
[1.168560, 0.148786, -0.054492, ..., -0.030141, -0.030610, 0.757597]

请注意，我们有几个向量，但最后一个是 Matrix 。我们必须转换smSplineEx到 python 分为两部分。

In [63]:

ro.r["smSplineEx1"].names
Out[63]:
<StrVector - Python:0x108525dd0 / R:0x1042ca7c0>
['time', 'y', 'y.true', 'all', 'Zt']
In [64]:

print com.convert_robj(ro.r["smSplineEx1"].rx(ro.IntVector(range(1, 5)))).head()
   time         y    y.true  all
1     1  5.797149  4.235263    1
2     2  5.469222  4.461302    1
3     3  4.567237  4.678477    1
4     4  3.645763  4.887137    1
5     5  5.094126  5.087615    1
In [65]:

print com.convert_robj(ro.r["smSplineEx1"].rx2('Zt')).head(2)
         0         1         2         3         4         5         6   \
1  1.168560  2.071261  2.944953  3.782848  4.584037  5.348937  6.078121   
2  0.148786  1.072013  1.948857  2.789264  3.593423  4.361817  5.095016   

         7         8         9     ...           88        89        90  \
1  6.772184  7.431719  8.057321    ...     0.933947  0.769591  0.619420   
2  5.793601  6.458153  7.089255    ...     0.904395  0.745337  0.599976   

         91       92        93        94        95        96        97  
1  0.484029  0.36401  0.259959  0.172468  0.102133  0.049547  0.015305  
2  0.468893  0.35267  0.251890  0.167135  0.098986  0.048026  0.014836  

[2 rows x 98 columns]

com.convert_robj(ro.r["smSplineEx1"])由于混合数据结构问题，将无法工作。