title: time series dimensionality reduction

text: I have a call center data (such as one below) that has call data collected every 15 minutes. For a day the periodicity is 96 and for a week the periodicity is (7 x 96 = 672). If I would like to forecast this data I can use models that handle multiple seasonality such as `DSHW, TBATS and BATS` in R forecast package. As you can imagine having a periodicity of 96 and 672 would require tremendous amount of computational resource and not very many techniques can handle such a long periodicity. I recently came across an excellent article by Shen that uses Singular Value Decomposition to reduce the dimension/periodicity of the time series, so we could have a parsimonious model and the model would be more efficient/require less computational resource. I'm open to other methods such as discrete wavelet transformation (dwt), discrete Fourier transformation (dft) etc.,                   Time in mins                         Date      0:15 0:20 0:45 1:00 1:15 1:30     1/1/2012    41  12  33  42  38  27     2/1/2012    50  12  20  14  27  46     3/1/2012    30  17  21  44  47  45     4/1/2012    26  15  10  13  50  38     5/1/2012    39  17  22  28  48  12     6/1/2012    30  10  25  30  21  21     7/1/2012    47  46  18  10  36  37     8/1/2012    26  37  41  36  37  12      **Below are my questions:**   1. How to implement SVD/DWT/DFT based dimensional reduction for the aforementioned time series data using tools like R and matlab ?   2. Are there any automatic procedures available for time series dimensional reduction ?


title: Regression factors and covariance matrix

text: I am trying to follow someone else's notes. They have two matrices. One is called `comfact` (company factors). This is a 580 x 5 matrix. The 580 rows represent 580 different companies. The 5 columns represent 5 betas for the 5 different factors that were used to model a companies growth rate. These beta's were calculated by OLS regression. The second matrix is called `covmatrix`. It is a 5 x 5 covariance matrix of the 5 factors above. There is line in the notes (below) which I'm not sure what it is doing and what the result is showing me?               comfact * covmatrix * comfact'


title: Interpretation of age coefficient in relative survival analysis

text: Suppose we have a Cox proportional hazards model with covariates including the patient's age, where the log hazard ratio for age is positive and significant, indicating increased risk for outcome in older patients. Then, with the same data and covariates we run a relative survival model. The (raw) coefficient for age now becomes negative and significant. The signs of the coefficients for the other covariates do not change. In this document, at the end of section 7.2, the negative sign is explained with the following comment: > It is not surprising that age is negatively associated with the outcome in > such analysis as older patients will be losing relatively less than the > younger ones. I don't understand this. Could someone expand this explanation to help me make more sense of it? I understand that relative survival models use the risk of death in the general population (using life tables) so that the relative survival ratio is the ratio of observed survival to expected survival in the general population. The example referred to in the above document, along with the Cox model can be obtained in R using               library(relsurv)     data(slopop)     data(rdata)          # cox model     summary(coxmodel <- coxph(Surv(time,cens)~age+sex+year, data=rdata))          # relative survival model     rstrans(Surv(time,cens)~age+sex+year+         ratetable(age=age*365.24,sex=sex,year=year),          data=rdata, ratetable=slopop)      which produces the following for the cox model:               coxph(formula = Surv(time, cens) ~ age + sex + year, data = rdata)            n= 1040, number of events= 547                      coef  exp(coef)   se(coef)      z Pr(>|z|)         age   5.938e-02  1.061e+00  4.324e-03 13.734  < 2e-16 ***     sex   9.011e-02  1.094e+00  9.585e-02  0.940  0.34716         year -2.736e-04  9.997e-01  9.177e-05 -2.981  0.00287 **       and this for the the relative survival model:               rstrans(formula = Surv(time, cens) ~ age + sex + year + ratetable(age = age *          365.24, sex = sex, year = year), data = rdata, ratetable = slopop)                    coef exp(coef) se(coef)     z       p     age  -0.013904     0.986 4.92e-03 -2.83 4.7e-03     sex   0.528578     1.697 1.01e-01  5.24 1.6e-07     year -0.000228     1.000 9.04e-05 -2.52 1.2e-02


title: cv.glmnet lambda stability

text: **Intro** I am running `cv.glmnet` from the `glmnet` package in `R`. I am running 10-fold cross validation 100 times on a dataset that has 25,000 observations and 150 variables. I am using the `parallel = TRUE` option, as that greatly speeds up results. **Problem Description**   * I notice that the minimum lambda value (`lambda.min`) and the lambda value within 1 standard error of the minimum (`lambda.1se`) are exactly the same in 99.9% of all 100 cross validation runs.    * Each cross validation run leads to the calculation of exactly 25 lambdas - no more and no less (recall that the default is that the number of lambdas depends on whether the deviance changes sufficiently from one run to the next). I have left the default settings as they are, which allows for the calculation of 100 lambdas unless the deviance stops changing significantly.   * The cross validation error varies minimally (average of 0.002).    * If I run another few iterations of 100 10-fold cross validations on a different split into training and test populations, I get a slightly different minimum lambda and lambda within 1 standard error of the minimum from the other iterations. The above points hold, however, for the given iteration (i.e. same lambda, same exact number of lambdas calculated per cross validation run, minimal change in cross validation error). **Final Question** Why are all my cross validation results within a given iteration so similar (perhaps because all folds in the cross validation are always very similar)?


title: Ascertaining what clustering algorithm to use in various situations

text: It is said that kMeans clustering works as long as we don't have clusters of differing   1. sizes,   2. densities,   3. and non-spherical shapes I understand how one might check the sizes and densities of data points, but I am unsure how one would check for non-spherical clusters, especially if one merely has a conceptual understanding of a dataset. So, **Q1.** Since most people with real world data sets would not know if they have non-spherical clusters, would it not be a mistake to apply kMeans? **Q2.** What would be a more suitable choice of a clustering algorithm if we only want to _partition_ our data set, given we don't know the shape of our data? _By partition, I refer to partitioning cluster algorithms such as kMeans and kMedoids._ **Q3.** What mind map would experienced clustering users have when they intend to apply partitioning cluster algorithms? For example, the following is mine, but I am inexperienced in clustering:   1. check if the data is comprised of continuous, or continous + discrete variables   2. pick euclidean/manhattan/mahalanobis distance measure for former, or gower for latter   3. perhaps standardise first and then apply the distance measure/or correlation if a similarity matrix is sought   4. apply a partitioning cluster algorithm; kMeans if we know our data shape and outliers are not a problem. kMedoids if we don't want to be affected by outliers ( **although** the definition of a centroid/prototype as found by kMeans seems to change if we apply kMedoids, i.e. a centroid in kMeans is understandably a prototype of each cluster, while in kMedoids this centroid is simply another data point and not a prototype)   5. validate cluster results (using the Elbow test/Dunn index/Silhouette Coeff if focusing only on internal measures) and also check for clustering tendency Am I missing something important in the above? Many thanks, and please pardon the ignorance.


title: Categorical Predictors in multiple regression

text: Simple question: I have a dataset, in which all the multivariate x variables (x0, x1, x2, x3..) are continuous, and all the y variables are categorical (distributed (equally) between 1-20 categories). Can this type of data be processed via multiple linear regression? example:               data = [2.3, 4.5, 0.7, 2.1], [2]            [3.4, 2.1, 0.7, 1.2], [5]            [3.1, 1.2, 4.5, 4.1], [4]            [4.1, 4.4, 3.2, 6.2], [11]            ...      Thanks. :-)


title: Comparing residuals between OLS and non-OLS regressions

text: Suppose you want to estimate a linear model: ($n$ observations of the response, and $p+1$ predictors) $$\mathbb{E}(y_i) = \beta_0 + \sum_{j=1}^p \beta_j x_{ij}$$ One way to do this is through the OLS solution, i.e. choose the coefficients so that the sum of square errors is minimum: $$(\beta_0,\beta_1,\cdots,\beta_p)^T = \underset{\beta_0,\beta_1,\cdots,\beta_p}{\arg \min} \sum_{i=1}^{n} \left( y_i - \beta_0 - \sum_{j=1}^p \beta_j x_{ij} \right)^2 $$ Alternatively, you could use another loss function, like the sum of the absolute deviations, so that: $$(\beta_0,\beta_1,\cdots,\beta_p)^T = \underset{\beta_0,\beta_1,\cdots,\beta_p}{\arg \min} \sum_{i=1}^{n} \left| y_i - \beta_0 - \sum_{j=1}^p \beta_j x_{ij} \right| $$ Suppose you have found the parameters for the two models, and want to choose the model with the smallest value of the loss function. How can you compare the minimum values attained by the loss functions in general? (i.e. not just this specific case - we could also try other $L_p$ based loss functions) There seems to be a difference in the scale of the functions - one deals with squares while the other does not.


title: Validating a multivariate categorical model

text: I assume that my population is a sample of an unknown multivariate categorical distribution $\mathbf{X} = (X_1, X_2, \ldots, X_k)$. From this population, a sample $\mathbf{X^*}$ is available, I assume $\mathbf{X^*} \sim \mathbf{X}$. I am evaluating different methods to reconstruct $\mathbf{X}$ from such a sample. What would be a good way to assess the "quality" of the reconstruction method from the sample alone, for both reproduction of the data and generalization error? I am thinking about applying cross validation with, say, $B=500$ repeated subdivisions, and computing an error measure (such as the square root of the mean square error) to obtain a distribution of the error. My questions are:   * When computing the error measure, what should be used as benchmark -- only the test data, or both training and test data? Does it make sense to compute two error measures -- one against the training, and the other against the test data?   * How should I split the sample? $90/10$, $1-e^{-1}/e^{-1}$, $50:50$, ...?   * Would it be more appropriate to use out-of-bootstrap error (=sampling with replacement)?   * Which error measure (=loss function?) is best? Any hints on terminology or literature are also more than welcome.


title: Likelihood of errors based on a sample

text: I frequently do something like: load a bunch of data, and then scan some fraction of it randomly to verify that no errors occurred. The more data I verify, the greater my certainty that no errors occurred anywhere. I'm curious about how I can express this formally. My thought process is: Suppose there are $n$ records, of which $k$ contain errors. The probability of a randomly selected one containing an error is $k/n$ and conversely, the probability of selecting one without an error is $1-k/n$. If I model this as a Bernoulli process, I believe the probability of $t$ samples all having no error is $\left(1-k/n\right)^t$. So given some information about the data, I can get the probability of finding no errors. But what I really want is the reverse: given that I found no errors, what confidence do I have that there are less than $k_{act}$ errors in actuality? I think I could use Bayes' theorem here, but that would require me having some priors which seem difficult to estimate. Can I just set my formula equal to, say, $.95$ and solve, and then claim I've found a 95% confidence interval? By the way, I'm interested for curiosity's sake. I think a more realistic attempt would take into account that errors would tend to propagate, for example.


title: Correlation matrix

text: Is it good practice to consistently use a correlation matrix when the direction of the relationship between dependent and independent variables is unknown (that is, when we’re uncertain about an inverse or direct relationship between variables)? Or, should established research be used as the basis for determining the type of relationship between variables?


title: Different ways to write interaction terms in lm?

text: I have a question about which is the best way to specify an interaction in a regression model. Consider the following data:               d <- structure(list(r = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("r1","r2"),          class = "factor"), s = structure(c(1L, 1L, 1L, 1L, 1L,           2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L),          .Label = c("s1","s2"), class = "factor"), rs = structure(c(1L, 1L,          1L,1L, 1L,2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L),         .Label = c("r1s1","r1s2", "r2s1", "r2s2"), class = "factor"),           y = c(19.3788027518437, 23.832287726332, 26.2533235300492,          15.962906892112, 24.2873740664331, 28.5181676764727, 25.2757801195961,          25.3601044326474, 25.3066440027202, 24.3298865128677, 32.5684219007394,          31.0048406654209, 31.671238316086, 34.1933764518288, 36.8784821769123,          41.6691435168277, 40.4669714825801, 39.2664137501106, 39.4884849591932,          49.247505535468)), .Names = c("r","s", "rs", "y"),           row.names = c(NA, -20L), class = "data.frame")      Two equivalent ways to specify the model with interactions are:               lm0 <- lm(y ~ r*s, data=d)     lm1 <- lm(y ~ r + s + r:s, data=d)      My question is if I could specify the interaction considering a new variable (rs) with the same levels of interaction:               lm2 <- lm(y ~ r + s + rs, data=d)      What advantages/disadvantages have this approach? And why the results from these two approaches are different?               summary(lm1)          lm(formula = y ~ r + s + r:s, data = d, x = TRUE)                 coef.est coef.se     (Intercept) 21.94     1.46       rr2         11.32     2.07       ss2          3.82     2.07       rr2:ss2      4.95     2.92       ---     n = 20, k = 4     residual sd = 3.27, R-Squared = 0.87               summary(lm2)          lm(formula = y ~ r + s + rs, data = d, x = TRUE)                 coef.est coef.se     (Intercept) 21.94     1.46       rr2         11.32     2.07       ss2          8.76     2.07   # ss2 coef is different from lm1     rsr1s2      -4.95     2.92       ---     n = 20, k = 4     residual sd = 3.27, R-Squared = 0.87


title: How should a math-paranoid programmer learn stats?

text: (this is the opposite of "introduction to statistics for mathematicians") I'm a programmer who loves writing logic and even a libsvm vector embedding library to my name. Now I'm getting asked questions about:   * "exogenous factors"   * "sample randomness"   * distributions, confidence intervals, etc. There are plenty of self-taught programmers, so how should I self-teach statistics?


title: Set up a classification experiment via libsvm (MATLAB)

text: I would like to set up an experiment in MATLAB, to predict the class of a set of text instances in a two-class problem (e.g., the text talks/does not talk about `basketball`). By doing this, I would like to draw **learning curves** by varying the number of samples `M` and consequently understand which is the best `M` to be used to avoid high bias/high variance problems. Consequently, I would like to compute $J_{train}(\theta)$, $J_{test}(\theta)$ and $J_{cv}(\theta)$ (respectively: training error, test error and cross-validation error) on my experiment. I created a question last week, but I guess it was misplaced and maybe it is better to ask it here. I'll try to subdivide my question in different subquestions. **Context** I looked for some useful documentation on `libsvm`, so as to understand how to use its parameters and set up properly the experiment. Unfortunately, all I can find is a list of random parameters without any further explanation, e.g.:               -d degree : set degree in kernel function (default 3)     -g gamma : set gamma in kernel function (default 1/num_features)     -r coef0 : set coef0 in kernel function (default 0)      Thus, as an initial action, I went for a classical setting, without specifying any classifier parameter:               model = svmtrain(trainingSet.labels, featureVectors, '-b 1');     [predict_label, accuracy, prob_estimates] = svmpredict(testSet.labels, testFeatureVectors, model, '-b 1');      where the feature vectors are integer vectors of length $L \sim 4000$ (TF-IDF approach, each cell indicates the number of occurrences of the specific word in a text instance). **Smoothness of test error** I tested the training-testing pipeline specified above by changing the number of samples $M$ that is fed into the classifier in training-testing phase. Then, I extracted the test error:               testError = accuracy(2);      When the number of results in the training set $M$ is increasing, I would expect the test error to change monotonically. However, what I obtain is the following figure (x-axis: $M$, y-axis: training error $J_{test}(\theta)$): ![test error](http://i.stack.imgur.com/FtSaQ.png) Is this a "normal" behavior? As I mentioned above, I was expecting something smoother than this. Thus, I see here different options:   1. either I am using the wrong parameters to set up the classifier   2. or I misinterpreted the meaning of the output   3. or I built wrongly the train-test pipeline Is there a more appropriate way of setting up the experiment? **Training error** How to compute the error $J_{train}(\theta)$ on the training set? I rerun the prediction on the **training set** and then extracted the MSE:               [~, accuracy, ~] = svmpredict(trainingSet.labels, featureVectors, model, '-b 1');     trainError(p) = accuracy(2);      However, again, the training error (or what I suppose is the training error) has this shape (which, again, is not monotonic): ![training error](http://i.stack.imgur.com/EgDUk.png) **Cross-validation** How to use cross-validation and consequently compute the cross-validation error $J_{cv}(\theta)$? I am using:               [~,accuracy.~] = svmtrain(crossValidationSet.labels, featureVectors);     CVerror = accuracy(2);      but again, as above, I don't know whether this makes sense. **Overall result** Here, I am drawing cross validation error (red) and training error (blue). Not monotonic. ![cross validation error](http://i.stack.imgur.com/HQUuv.png)


title: Mediation with interaction and 3-level categorical IVs using R 'mediation' package

text: Using the ‘mediation’ package by Tingley and colleagues in R, I’m trying to conduct a bootstrapped mediation analysis for an interaction between a 2-level categorical IV (‘condstorm’) and another 3-level categorical IV (‘condmean’). The mediator (‘trustmed2’) is dichotomous, and the dependent variable (‘AttInd’) is continuous. Categorical 2-level IV ‘condstorm’               table(d$condstorm)              nostorm   storm          604     570       Categorical 3-level IV ‘condmean’               table(d$condmean)          both high mean      413  375  386       Dichotomous mediator ‘trustmed2’               table(d$trustmed2)              0   1          305 866      Continuous DV ‘AttInd’               summary(d$AttInd)     Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's      0.0000  0.3571  0.6200  0.5801  0.8214  1.0000       3       The code and output for the model with the mediator is:               contrasts(d$condstorm)         storm     nostorm     0     storm       1          contrasts(d$condmean)          nonevhigh nonevboth     both         0         1     high         1         0     mean         0         0               med.fit <- glm(trustmed2 ~ condmean*condstorm + gwbelief + trustmed1 + coastdwell +                        partyID + polor + gender + age + race + educ + income +                        region, weights=weight1, data = d, family = binomial('probit'))          summary(med.fit)          Call:     glm(formula = trustmed2 ~ condmean * condstorm + gwbelief + trustmed1 +          coastdwell + partyID + polor + gender + age + race + educ +          income + region, family = binomial("probit"), data = d, weights = weight1)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.9982  -0.1429   0.2196   0.4594   3.2721            Coefficients:                                       Estimate Std. Error z value Pr(>|z|)         (Intercept)                      -1.876115   0.284494  -6.595 4.27e-11 ***     condmeannonevhigh                 0.170809   0.175454   0.974  0.33029         condmeannonevboth                 0.497041   0.181383   2.740  0.00614 **      condstormstorm                    0.285569   0.182081   1.568  0.11680         gwbeliefnotvref                   0.056512   0.631022   0.090  0.92864         gwbeliefnotvprobs                 0.746556   0.129610   5.760 8.41e-09 ***     trustmed1                         1.908566   0.121064  15.765  < 2e-16 ***     coastdwellyes                     0.320369   0.128445   2.494  0.01262 *       partyIDvdem                       0.585477   0.148862   3.933 8.39e-05 ***     partyIDvref                      -0.468275   0.654078  -0.716  0.47403         partyIDvrep                       0.058146   0.130621   0.445  0.65621         polorvref                        -0.316949   0.508007  -0.624  0.53269         polorvlib                        -0.058482   0.160609  -0.364  0.71577         polorvcons                       -0.369565   0.122246  -3.023  0.00250 **      genderMale                       -0.009229   0.105253  -0.088  0.93013         agev1824                          0.235127   0.194856   1.207  0.22756         agev2534                          0.099723   0.178256   0.559  0.57587         agev3544                          0.119399   0.169987   0.702  0.48243         agev4554                          0.221921   0.176420   1.258  0.20842         agev5564                          0.173343   0.170355   1.018  0.30890         racevhisp                         0.033141   0.177162   0.187  0.85161         racevblack                       -0.323867   0.188305  -1.720  0.08545 .       racevother                        0.475179   0.234313   2.028  0.04256 *       educhsvsome                      -0.002255   0.132574  -0.017  0.98643         educhsvcollgrad                   0.201980   0.146240   1.381  0.16723         educhsvref                        1.304337   0.666427   1.957  0.05032 .       incomev3050                      -0.095963   0.159290  -0.602  0.54688         incomev5075                       0.184528   0.168946   1.092  0.27473         incomev75100                      0.028269   0.191474   0.148  0.88263         incomev100                       -0.249873   0.167648  -1.490  0.13610         regionvmw                        -0.058581   0.163619  -0.358  0.72032         regionvne                         0.315495   0.177227   1.780  0.07505 .       regionvs                          0.204790   0.145384   1.409  0.15895         condmeannonevhigh:condstormstorm -0.061962   0.263472  -0.235  0.81407         condmeannonevboth:condstormstorm -0.665467   0.253958  -2.620  0.00878 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 1316.98  on 1170  degrees of freedom     Residual deviance:  767.03  on 1136  degrees of freedom       (3 observations deleted due to missingness)     AIC: 836.49          Number of Fisher Scoring iterations: 6      And the code and output for the model with the mediator included as a predictor on the dependent variable:               out.fit <- lm(AttInd ~ trustmed2 + condmean*condstorm + trustmed1 + gwbelief + coastdwell +                       partyID + polor + gender + age + race + educ + income +                       region, weights=weight1, data = d)          summary(out.fit)          Call:     lm(formula = AttInd ~ trustmed2 + condmean * condstorm + trustmed1 +          gwbelief + coastdwell + partyID + polor + gender + age +          race + educ + income + region, data = d, weights = weight1)          Weighted Residuals:          Min       1Q   Median       3Q      Max      -0.79365 -0.12054  0.00739  0.12351  0.72377           Coefficients:                                        Estimate Std. Error t value Pr(>|t|)         (Intercept)                       1.754e-01  3.411e-02   5.143 3.18e-07 ***     trustmed2                         2.777e-01  1.940e-02  14.313  < 2e-16 ***     condmeannonevhigh                -8.754e-03  2.174e-02  -0.403 0.687311         condmeannonevboth                 2.548e-02  2.148e-02   1.186 0.235756         condstormstorm                    3.977e-02  2.250e-02   1.767 0.077433 .       trustmed1                         2.084e-02  2.014e-02   1.035 0.301053         gwbeliefnotvref                  -1.063e-01  9.071e-02  -1.172 0.241589         gwbeliefnotvprobs                 1.574e-01  1.818e-02   8.656  < 2e-16 ***     coastdwellyes                     1.100e-02  1.485e-02   0.740 0.459293         partyIDvdem                       6.283e-02  1.618e-02   3.884 0.000109 ***     partyIDvref                      -2.558e-01  8.716e-02  -2.934 0.003408 **      partyIDvrep                      -4.893e-02  1.688e-02  -2.898 0.003829 **      polorvref                         8.509e-02  6.617e-02   1.286 0.198751         polorvlib                         6.199e-02  1.742e-02   3.558 0.000389 ***     polorvcons                       -7.015e-03  1.571e-02  -0.447 0.655204         genderMale                       -5.778e-03  1.267e-02  -0.456 0.648363         agev1824                          4.624e-03  2.391e-02   0.193 0.846692         agev2534                         -1.773e-02  2.222e-02  -0.798 0.425121         agev3544                          2.952e-02  2.107e-02   1.401 0.161448         agev4554                          1.590e-02  2.159e-02   0.736 0.461613         agev5564                         -9.967e-03  2.073e-02  -0.481 0.630765         racevhisp                         4.936e-02  2.030e-02   2.432 0.015174 *       racevblack                       -1.585e-02  2.224e-02  -0.713 0.476104         racevother                        4.969e-02  2.518e-02   1.973 0.048687 *       educhsvsome                      -1.016e-02  1.612e-02  -0.630 0.528528         educhsvcollgrad                  -1.209e-02  1.713e-02  -0.706 0.480595         educhsvref                        2.383e-01  8.213e-02   2.901 0.003791 **      incomev3050                      -9.128e-04  1.964e-02  -0.046 0.962941         incomev5075                       1.522e-02  1.996e-02   0.762 0.446132         incomev75100                     -3.185e-02  2.280e-02  -1.397 0.162813         incomev100                        1.344e-02  2.013e-02   0.668 0.504422         regionvmw                         3.467e-02  2.027e-02   1.710 0.087541 .       regionvne                         1.861e-02  2.015e-02   0.924 0.355812         regionvs                         -7.531e-05  1.716e-02  -0.004 0.996499         condmeannonevhigh:condstormstorm -4.354e-03  3.152e-02  -0.138 0.890175         condmeannonevboth:condstormstorm -5.284e-02  3.090e-02  -1.710 0.087514 .       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.2103 on 1135 degrees of freedom       (3 observations deleted due to missingness)     Multiple R-squared:  0.4546,    Adjusted R-squared:  0.4378      F-statistic: 27.03 on 35 and 1135 DF,  p-value: < 2.2e-16      The model without the mediator is:               out.fit <- lm(AttInd ~ condmean*condstorm + trustmed1 + gwbelief + coastdwell +                      partyID + polor + gender + age + race + educ + income +                       region, weights=weight1, data = d)          summary(out.fit)          Call:     lm(formula = AttInd ~ condmean * condstorm + trustmed1 + gwbelief +          coastdwell + partyID + polor + gender + age + race + educ +          income + region, data = d, weights = weight1)          Weighted Residuals:         Min      1Q  Median      3Q     Max      -1.1727 -0.1371  0.0130  0.1366  0.8077           Coefficients:                                        Estimate Std. Error t value Pr(>|t|)         (Intercept)                       0.1940313  0.0370161   5.242 1.89e-07 ***     condmeannonevhigh                -0.0019461  0.0236085  -0.082  0.93432         condmeannonevboth                 0.0474046  0.0232707   2.037  0.04187 *       condstormstorm                    0.0527524  0.0244159   2.161  0.03094 *       trustmed1                         0.1833521  0.0180674  10.148  < 2e-16 ***     gwbeliefnotvref                  -0.1093369  0.0985135  -1.110  0.26729         gwbeliefnotvprobs                 0.2074035  0.0193781  10.703  < 2e-16 ***     coastdwellyes                     0.0241253  0.0161011   1.498  0.13432         partyIDvdem                       0.0858665  0.0174835   4.911 1.04e-06 ***     partyIDvref                      -0.2999944  0.0945990  -3.171  0.00156 **      partyIDvrep                      -0.0468268  0.0183366  -2.554  0.01079 *       polorvref                         0.0560635  0.0718304   0.780  0.43526         polorvlib                         0.0568556  0.0189179   3.005  0.00271 **      polorvcons                       -0.0293808  0.0169730  -1.731  0.08372 .       genderMale                       -0.0072970  0.0137554  -0.530  0.59588         agev1824                          0.0162995  0.0259507   0.628  0.53007         agev2534                         -0.0134314  0.0241301  -0.557  0.57790         agev3544                          0.0326065  0.0228762   1.425  0.15433         agev4554                          0.0281826  0.0234333   1.203  0.22935         agev5564                         -0.0042263  0.0225102  -0.188  0.85111         racevhisp                         0.0484325  0.0220418   2.197  0.02820 *       racevblack                       -0.0292002  0.0241349  -1.210  0.22658         racevother                        0.0682464  0.0273067   2.499  0.01259 *       educhsvsome                      -0.0130208  0.0175088  -0.744  0.45723         educhsvcollgrad                  -0.0052682  0.0185980  -0.283  0.77703         educhsvref                        0.3522019  0.0887800   3.967 7.73e-05 ***     incomev3050                      -0.0065261  0.0213276  -0.306  0.75967         incomev5075                       0.0200734  0.0216792   0.926  0.35468         incomev75100                     -0.0336651  0.0247648  -1.359  0.17429         incomev100                       -0.0006802  0.0218303  -0.031  0.97515         regionvmw                         0.0293874  0.0220138   1.335  0.18216         regionvne                         0.0317515  0.0218553   1.453  0.14655         regionvs                          0.0085483  0.0186236   0.459  0.64632         condmeannonevhigh:condstormstorm -0.0068725  0.0342330  -0.201  0.84092         condmeannonevboth:condstormstorm -0.0860072  0.0334623  -2.570  0.01029 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.2284 on 1136 degrees of freedom       (3 observations deleted due to missingness)     Multiple R-squared:  0.3562,    Adjusted R-squared:  0.3369      F-statistic: 18.49 on 34 and 1136 DF,  p-value: < 2.2e-16      We want to examine mediation in this interaction effect. We have tried using the code:               mediate(med.fit, out.fit,  treat = 'condstorm', control.value = 'nostorm', treat.value = 'storm', mediator = 'trustmed2', covariates=list(condmean='both'), boot=TRUE, sims=5000)      Which gave the following output:               Causal Mediation Analysis           Confidence Intervals Based on Nonparametric Bootstrap          (Inference Conditional on the Covariate Values Specified in `covariates')          Estimate 95% CI Lower 95% CI Upper p-value          ACME           -0.017705    -0.036626    -0.003284    0.04     ADE            -0.014082    -0.045811     0.023155    0.53     Total Effect   -0.031788    -0.069698     0.012371    0.14     Prop. Mediated  0.516591     1.174065     8.857665    0.04          Sample Size Used: 1171           Simulations: 5000      But this mediation only compares respondents at the level of “both” for the condmean variable and excludes other conditions, so examines the effect of "Both Condmean x Condstorm" MINUS the effect "condstorm", which is -.08 - (.05) = -.03 (the total effect in the output above). The total effect should instead be the same as the coefficient of "condmeannonevboth:condstormstorm" in the second interaction term in the model, which is -.086, p < .01. Is it possible for us to examine this effect using the ‘mediation’ package in R? Any help would be much appreciated!


title: Why is this distribution uniform?

text: We are investigating Bayesian statistical testing, and come across an odd (to me atleast) phenomenon. Consider the following case: we are interested in measuring which population, A or B, has a higher conversion rate. For a sanity check, we set $p_A = p_B$, that is, the probability of conversion is equal in both groups. We generate artificial data using a binomial model, e.g. $$n_A \sim \text{Binomial}(N, p_A)$$ We then try to estimate the $p_A, p_B$ using a Bayesian beta-binomial model so we get posteriors for each conversion rate, e.g. $$P_A \sim \text{Beta}(1 + n_A, N - n_A +1 )$$ Our test statistic is computed by calculating $S = P(P_A > P_B\; |\; N, n_A, n_B)$ via monte carlo. What surprised me was that if $p_A = p_B$, then $S \sim \text{Uniform(0,1)}$. My thoughts were that it would be centered around 0.5, and even converge to 0.5 as the sample size, $N$, grows. My question is, why is $S \sim \text{Uniform(0,1)}$ when $p_A = p_B$? * * * Here's some Python code to demonstrate:               %pylab     from scipy.stats import beta     import numpy as np     import pylab as P          a = b = 0.5     N = 10000     samples = [] #collects the values of S     for i in range(5000):         assert a==b         A = np.random.binomial(N, a); B = np.random.binomial(N, b)         S = (beta.rvs(A+1, N-A+1, size=15000) > beta.rvs(B+1, N-B+1, size=15000)).mean()          samples.append(S)          P.hist(samples)     P.show()


title: How to compare pairs of coefficients within a glmm with binomial error

text: I have a generalised linear mixed model with 34 explanatory variables (over 130,000 observations for each). 10 of these variables are different unprotected habitat types, and another 10 are the same habitat types but protected for conservation. (The factor for the random effect E.V. is region). I want to compare whether the coefficients (within this single model) for each habitat type when protected or not protected are significantly different i.e. compare beta(broadleaf.protected) with beta(broadleaf.unprotected), beta(conifer.protected) with beta(conifer.unprotected) etc. to find out whether they have a significantly different effect on my D.V. (visit probability). My explanatory variables are standardised to mean 0 sd 1.


title: Effect of using regressor with low variance in a fixed effects regression

text: Does including a regressor with minimal variance behave analogously to including an additional fixed effect in the regression? If yes, what could be the implications for that? I am using a panel dataset of 29 cross-sectional units in 47 years.


title: Dealing with missing data - glmer in lme4 package

text: Dear StackExchange community, I have an unbalanced data set / data set with missing values, consisting of 20 submersible acoustic receivers that have been range tested on 8 days (Both receiver ID and Day are treated as random effects in my model). My aim is to test the effects of multiple environmental variables on the detection range of these receivers. Unfortunately, due to technical issues, a total of 5 receivers could not be tested every day (e.g., 2 receivers could not be tested on 2 days, 3 could not be tested on 1 day). What would be the right way to proceed? I have 2 choices: 1) reduce data set (exclude all receivers that are not tested on every day), or 2) work with full data set, including missing values. The first choice is easy, but may not be the best choice (as far as I've read online). besides, collected data will be thrown away.. The second choice however seems more difficult to work with. I read that glmm in the lme4 package can deal with missing values, however, the only thing it does is automatically exclude all rows that contain NAs. So let's say I choose the second option, and let the model run with missing values that automatically get deleted. How would this affect hypothesis testing? In other words, is the interpretation of p-values just as straightforward as when it would be a balanced design? [EDIT: I worked out the full data analysis that excluded those receivers for my masters project in order to maintain a balanced data set, however for a publication I would like to analyse my data using the second choice, as I think it's a better one. There's not much literature surrounding this subject as far as I'm aware, hence my post to this forum.] Thank you for your consideration, Maurits


title: Least squares with multiple constraints

text: I work with a regression with ARMA errors and I want to use LASSO to shrink the coefficients and select my variables. This topic is discussed in the article Wu et al. (2012). So, the problem I have to minimize is: $$ L(\xi)=||Y-X\xi||^2+\lambda_1\sum_{i=1}^{d}\frac{|\beta_i|}{|\beta_{i}^{OLS}|}+\lambda_2\sum_{i=1}^{p}\frac{|\phi_i|}{|\phi_{i}^{OLS}|}+\lambda_3\sum_{i=1}^{q}\frac{|\theta_i|}{|\theta_{i}^{OLS}|} $$ where: $\xi=[\beta,\phi,\theta]$ and $\xi^{OLS}=[\beta^{OLS},\phi^{OLS},\theta^{OLS}]$ is the least squares estimate. Do you know if there is a function in R that can solve this minimization? Something similar with glmnet / l1ce / lars, but with multiple constraints on the parameters? Thank you!


title: Homogeneity testing of baseline characteristics in medical trials

text: I have been reading medical journals and they repeatedly show baseline characteristics of samples from a randomised controlled trial, which they have then tested to ensure no differences between the two groups under study. For example in the operation group you have 39 males and 3 females. And in the non-op group you have 35 males and 8 females. The males would have a p-value, and then the females would have a p-value. I was wandering what exactly is being tested and how. My initial thoughts were a chi-squared test, although this would produce only 1 p-value, and this would be looking at the distribution of patients between the two groups. The papers seem to suggest they assess just the females distribution, and then separately the males distribution. But is this right and how would one do it?


title: GLMM - time or day/date as a random factor

text: I have searched for many days trying to find the answer to this question, and am still not 100% sure I am happy with my conclusion. I am interested in looking at the effects of environmental variables on detectability of marine mammals. I have Animal Counts as my response variable, and have Cloud cover, sea state, glare level, etc. as environmental variables. These were measured throughout the day over good weather days in a season. This was done only at one location. I am not interested in the effect of time of day or Date (day in the season) on counts, although these will have an influence (and will counts over these variables are likely to be auto-correlated). Rather, I am interested in the effect of environmental conditions on the numbers counted. My thought was to run a GLMM with Time in the day (maybe counts over hourly Time Blocks) and Date as random factors (Time Block nested under Date), and then include the environmental variables as fixed factors. So Date would become the 'subject' in the GLMM. I realise that opinion on what should be considered as a random variable varies, and there are a range of approaches to a single problem. I have seen a paper where it appears that Day itself was used as a subject. Any thoughts on whether this approach should be appropriate or suggested better approaches are much appreciated. Hopefully this will be useful to others as well. Chandra


title: Probability of getting ace after getting a king

text: Consider a regular 52 card deck of playing cars. Suppose 2 cards are drawn at random from the deck with-out replacement. What is the probability that the second card is an ace, given that the first card is a king?


title: Computing c-index for an external validation of a Cox PH model with R

text: First off, I'll state that I'm aware many questions get asked about the c-index. I've searched this site and others, and I haven't found an answer for my situation. I can successfully use `validate()` in the `rms` package to calculate the Dxy and c-index for my boot-strapped internal validation. Now I need a c-index for my external validation. I just externally validated my model with an independent data set, and produced the desired plot of actual vs. predicted probabilities using `val.surv()`. Unfortunately the abstract I'm submitting cannot have figures, so I feel compelled to report a c-index for the external validation. I have searched this site and the R help archive and haven't found a conclusive answer for how to calculate the c-index on external validation. I've seen mentions of using `rcorr.cens()` in the `Hmisc` package, but it seems to me that you can only use it for a single variable's concordance, not an entire model. As of yet I can't find a way to use `val.surv()` for calculating a c-index as well. I posted some sample code below, including a test data set to act like an external validation set. I really appreciate your help in calculating a c-index from an external validation with an independent data set.               library(rms)     library(Hmisc)     data(veteran)          ##Create a Cox PH model for the training data.     survmod=with(veteran,Surv(time,status))     cox.mod=cph(survmod~celltype+karno,data=veteran,x=T,y=T,surv=TRUE,time.inc=5*365)          ##Here is the test data set that is the external "independent" data.     test_dat=data.frame(trt=replicate(500,NA), celltype=replicate(500,NA), time=replicate(500,NA), status=replicate(500,NA), karno=replicate(500,NA), diagtime=replicate(500,NA), age=replicate(500,NA), prior=replicate(500,NA))     for(i in seq(8)){     test_dat[,i]=sample(veteran[,i],500,replace=T)     }          ##Validate the model with the test data     test_surv=with(test_dat,Surv(time,status))     validated=val.surv(cox.mod,newdata=test_dat,S=test_surv)          ##Now what I need is to take the external validation and compute the      #c-index.  This is where I'm stuck.       #I've seen people mention `rcorr.cens()`, but I can't figure out a way to use      #`rcorr.cens()` with a Cox model of several variables.  I appreciate your help!


title: GEEs in case of a small number of clusters with strongly heterogeneous cluster size

text: My data set includes 400 records. Each record comprises values for the binary outcome variable $y$ and 12 categorical predictor variables $x_1, ..., x_{12}$, most of which are binary too. The records originate from 10 different studies, of which one is much larger than the remaining ones (it contributes nearly three-third of the records, while each of the remaining studies contributes 4 to 30 records). My plan is to fit a logistic regression model. Since I cannot rule out that records from the same study are correlated, I consider using GEEs instead of ordinary logistic regression modeling. The records from the same study would form the clusters. As working correlation structure I would use 'exchangable'. However, I am still unsure if GEEs are a suitable approach in my case, which is distinguished by   1. the low number of clusters (viz., ten)    2. the large differences in cluster size.  Are there any advices/experiences on this? I found in the literature that actually 40, 50 or even more clusters would be needed (opinions seem to diverge) in order to get reliable 'sandwich' errors for the model parameter estimates. For smaller numbers of clusters, some correction would be needed. However, I neither found much on how to precisely apply this correction, nor on the second issue, if differences in cluster size affects the results of GEEs. I would highly appreciate any advice.


title: Goodness of fit – Testing Cox proportional hazard assumption in R

text: I am performing survival analysis on credit data. I created a simple model with using interest rate:   `cox <\- coxph(Surv(periods,charged_off) ~ int_rate, data=notes)` I assumed that `int_rate` was a time-independent variable, but the following test rejects HA:               > cox.zph(cox)                 rho chisq        p     int_rate 0.0446  14.2 0.000169      **Same result for other variables such as loan amount:**               > cox <- coxph(Surv(periods,charged_off) ~ int_rate + loan_amnt, data=notes)     > cox.zph(cox)                  rho chisq        p     int_rate  0.0364  9.31 2.28e-03     loan_amnt 0.0317  8.84 2.95e-03     GLOBAL        NA 26.28 1.97e-06      **Plot for`int_rate`:** ![enter image description here](http://i.stack.imgur.com/Bz1oC.png) Why would these covariates be considered time dependent? Am I doing something wrong? Thanks.


title: Time Series Ahead Prediction in Neural Network, Large Scale Iterative Training

text: I am having trouble in implementing neural network to predict N points ahead. My only feature is previous time. I used elman recurrent neural network and also newff. In my scenario I need to predict 90 points ahead. First how I separated my input and target data manually: For Example:               data_in = [1,2,3,4,5,6,7,8,9,10]; //imagine 1:10 only defines the array index values.     N = 90; %predicted second ahead.     P(:, :)         T(:)     it could also be(2 theta time)  P(:, :)         T(:)      [1,2,3,4,5]    [5+N]              |                     [1,3,5,7,9]     [9+N]     [2,3,4,5,6]    [6+N]              |                     [2,4,6,8,10]    [10+N]     ...      until it reaches to end of the data > I have 100 input points and 90 output points in Elman recurrent neural > networks. What could be the most efficient hidden node size?               input_layer_size = 90;       NodeNum1 =90;           net = newelm(threshold,[NodeNum1 ,prediction_ahead],{'tansig', 'purelin'});     net.trainParam.lr       = 0.1;      net.trainParam.goal     = 1e-3;       //At the beginning of my training I filter it with kalman, normalization into range of [0,1] and after that I shuffled the data. 1) I won't able to train my complete data. First I tried to train complete M data which is around 900,000, which didn't gave me a solution. 2) Secondly I tried iteratively training. But in each iteration the new added data is merged with already trained data. After 20,000 trained data the accuracy start to decreases. First trained 1000 data perfectly fits in training. But after when I start iterativelt merge the new data and continue to training, the training accuracy drops very rapidly 90 to 20. For example.               P = P_test(1:1000) T = T_test(1:1000) counter = 1;          while(1)           net = train(net,P,T, [], [] );%until it reaches to minimum error I train it.           [normTrainOutput]    = sim(net,P,         [], [] );                P = [ P P(counter*1000:counter*2000)]%iteratively new training portion of the data added.          counter = counter + 1; end      This approach is very slow and after a point it won't give any good resuts. My third approach was iteratively training; It was similar to previous training but in each iteration, I do only train the 1000 portion of the data, without do any merging with previous trained data.For example when I train first 1000 data until it gets to minimum error which has >95% accuracy. After it has been trained, when I have done the same for the second 1000 portion of the data;it overwrites the weight and the predictor mainly behave as the latest train portion of the data.               > P = P_test(1:1000) T = T_test(1:1000) counter = 1;            while(1)     >       net            = train(net,P,T, [], [] ); % I did also use adapt()     >       [normTrainOutput]    = sim(net,P,         [], [] );     >         >       P = [ P(counter*1000:counter*2000)]%iteratively only 1000 portion of the data is added.     >     counter = counter + 1;      end      I wasn't able to successes any of this approaches. In each iteration I also observe that slight change on the alpha completely overwrites to already trained data and more focus onto the currently trained data portion. I won't able to come up with a solution to this problem. In iterative training should I keep the learning rate small and number of epochs as small. And I couldn't find an efficient way to predict 90 points ahead in time series. Any suggestions that what should I do to do in order to predict N points ahead, any tutorial or link for information. What is the best way for iterative training? On my second approach when I reach 15 000 of trained data, training size starts suddenly to drop. Iteratively should I change the alpha on run time? Any suggestion or the things I am doing wrong would be very appreciated. I also implemented recurrent neural network. But on training for large data I have faced with the same problems.


title: How to test for differences between groups over time when N is less than the number of levels

text: I have done a live cell imaging time course over 24 hours, and have a result for each hour. I have 3 experimental groups and 1 control group. What I want to know is if any of the experimental groups are significantly different from the control group, and if so, at what time does this first occur? I have 3 experimental replicates. I thought what I needed to do was a mixed design two way ANOVA, but it turns out (in SPSS at least) that I can't test for sphericity or homogeneity of variance if my N is less than my number of levels. How do I deal with this situation?


title: Model selection with nonlinear fitting? Statistical tests seem ambiguous

text: I'm working on fitting an exponential model $\mathrm{Flux} = A+Bt+F\left(\exp(t_0-t/T_r) + \exp(t-t_0/T_f)\right)^{-1}+...$ to astronomical data (a light curve). $A$, $B$, $F$, $t_0$, $T_r$, and $T_f$ are all free parameters in the model. I am trying to determine if a two-peak, three-peak, or four-peak model is appropriate for my dataset (with support from statistical tests). For each n-peak model (n=2,3,4) there is a corresponding exponential function in my model with it's own values of $t_0$ and F. I'm using nls in R to perform the fit. The problem I'm having is that it is ambiguous whether a 3-peak or 4-peak model is preferential. I also looked at a 5 peak model just for curiosity's sake. So far I have looked at a likelihood ratio test, the AIC and BIC and none seem to put preference on one model versus the other. Likelihood ratio test   _Using the 2-peak model as my null hypothesis_   3 peak: p=1.325564e-10   4 peak: p=2.220446e-16   5 peak: p=1.110223e-16   Here I see it as the 4-peak model being preferred over the 2-peak model, as the $H_0$ is rejected with a greater significance. But when you get p values this small, do differences this small really matter? _3-peak model as $H_0$_ 4 peak: 5.164630e-08   5 peak: 1.877226e-08   My interpretation: 4-peak model preferred. _4-peak model as $H_0$_   5 peak: .01661512   My interpretation: 5-peak model preferred over 4-peak model. AIC   2 peak:-58011.87   3 peak:-58053.36   4 peak:-58082.92   5 peak:-58087.11   My interpretation: 4-peak and 5-peak models equally preferred. But I would select the 4-peak model over the 5-peak model because it has fewer free parameters. But I question this because the difference between the AIC values (even between the 2-peak and 5-peak models) is so small. BIC   2 peak: -57988.89   3 peak: -58025.27   4 peak: -58049.72   5 peak: -58048.80   My interpretation: 4-peak model and 5-peak model equally preferred. Again, I would select the 4-peak model over the 5-peak model because it has fewer free parameters (also the BIC value is slightly smaller) I understand that when performing a LRT on nonlinear models it is preferable to use a non-parametric fit for your null hypothesis, however I have been unable to get that to work in R. I think my dataset may be too small (I have 95 data points). So my question really boils down to: have I interpreted these statistical tests properly. Also, are there any methods better than the ones I've already utilized to select one model over the other? Just for some context, the reason the model selection is so vital is that the values for $T_f$ under the 3-peak versus 4-peak model are very different and point to two differing physical explanations.


title: Unequal sample size t-test, but with group-effect

text: I am researching how people watch video lectures together. I have two conditions as follows: (1) 3 groups of 4 learners watch one video lecture on a shared display, but there is only one control (to pause for jump in video) (2) 3 groups of 4 learners watch one video lecture separately on individual displays, with individual control. They watch the same video on their own pace, but they can of course talk during watching, so there is a group effect. In the end, I want have a research question which is formulated by the following null hypothesis: **The two conditions have no effect on the frequency of pauses made.** This is a typical question that can usually be answered by doing a t-test (assuming normality and homogeneity). But the problem is that, for those who watch with single control and display, a group has a single log of interactions, however, for those who watch with multiple control and displays, a group has N (number of users) logs of interactions. So the sample size is different, but some samples are not independent, so we cannot directly use Mann-Whitney Test to replace t-test. I mean, the log files become something like this:                 group    user    numOfPauses  condition         A       A          20          Shared         B       A          32          Shared         C       A          22          Shared         D       D1         10          Individual         D       D2         23          Individual         D       D3         15          Individual         D       D4         3           Individual         E       E1         16          Individual         E       E2         38          Individual         E       E3         11          Individual         E       E4         9           Individual         F       F1         32          Individual         F       F2         21          Individual         F       F3         12          Individual         F       F4         17          Individual      Finally the number of observations from individually watching groups are N times the shared watching groups. Is there a way to balance the unequal size considering the group-effect? Or it is better to just make an average for individual watching groups, or expanding shared groups by duplicating the same value N times?


title: Markov chain model likelihood ratio test

text: Suppose I am using two Markov Chain Models, one with order $k=1$ and a second one with order $k=2$. I am "reducing" the higher order model to a $k=1$ model in order to have easier calculation possibilities. I train each model on the same data and also calculate the log likelihoods on the same data. Now I want to determine the log likelihood ratio test in order to make a model selection, as they are nested. To do so I need the LRT (which is straight forward) and the degrees of freedom. Currently, I am determining the df by calculating the difference between the parameters of the $k=2$ ($m^2(m-1)$) and the null model $k=1$ ($m(m-1)$). The problem now is that the degrees of freedom are very, very high and so I come up with a high p value all the time, which says that I should stick with my null model. I am unsure, if this is the right way to do so. The second order model is much sparser, so do I really need to calculate the worst case number of parameters, or can I make any limitations to that? Maybe someone can help me out with that. Cheers!


title: Data Exploration - How To Bin Data?

text: There are 1100 data points split into about 35 discrete groups and distributed towards the left. After running a logistic regression which failed, then running it again by only predicting based on subsets of the data which worked, it seemed likely that there is a nonlinear relationship. I am trying to determine the nature of this relationship to decide on how best to build my predictive model. The first image below is each discrete value on the x axis, and the y axis is the percent of target variable success. Since the first image doesn't reveal much I tried binning the data, so the three graphs below are the x variables binned into different groupings on the x axis, and the y axis is percent of the target variable success. It seems pretty hard to draw any conclusions as all these graphs show different things, or they don't show anything obvious, at least to me. Question 1- How do I best decide on how to bin data, to show relationships? Question 2- Are there any important relationships based on the graphs below that you can see? ![enter image description here](http://i.stack.imgur.com/NXbDe.jpg) ![enter image description here](http://i.stack.imgur.com/towPm.jpg)


title: Computing posterior distribution of bayesian lasso

text: I have a model: $$ \mathbf{y} = \mu\mathbf{1}_n + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon} $$ where $\boldsymbol{\epsilon} \sim N(\mathbf{0},\sigma^2\mathbf{I}_n)$. I have a joint prior: $$ \pi(\boldsymbol{\beta}, \sigma^2, \mu) = \pi(\mu) \pi(\sigma^2) \prod\limits_{j=1}^{p}\frac{\lambda}{2\sqrt{\sigma^2}}e^{-\lambda|\beta_j|/\sqrt{\sigma^2}} $$ where $\pi(\mu) \propto 1$. I want to compute joint posterior $\pi(\boldsymbol{\beta}, \sigma^2, \mu |\mathbf{y})$ and then marginalize out $\mu$. MY SOLUTION: According to me, likelihood function is: $$f(\mathbf{y} |\boldsymbol{\beta}, \sigma^2, \mu) = \frac{1}{(2\pi\sigma^2)^{\frac{n}{2}}} \exp\left( -\frac{1}{2\sigma^2}\parallel \mathbf{y}- \mu \mathbf{1}_n -\mathbf{X}\boldsymbol{\beta}\parallel_2^2 \right)$$ Using the Bayes theorem, I got the conclusion that joint log-posterior is proportional to: $$ \underbrace{\ln[\pi(\mu)]}_{\to \text{const.}} + \ln[\pi(\sigma^2)]-\frac{p+n}{2}\ln(\sigma^2) - \lambda \parallel \beta \parallel_1 \frac{1}{\sqrt{\sigma^2}} - \frac{1}{2\sigma^2}\parallel \mathbf{y}- \mu\mathbf{1}_n - \mathbf{X}\boldsymbol{\beta} \parallel_2^2 $$ HOWEVER, in paper about Bayessian Lasso (2008) written by Park and Casella they got: $$ \ln[\pi(\sigma^2)]-\frac{p+n-1}{2}\ln(\sigma^2) - \lambda \parallel \beta \parallel_1 \frac{1}{\sqrt{\sigma^2}} - \frac{1}{2\sigma^2}\parallel \mathbf{\tilde{y}}-\mathbf{X}\boldsymbol{\beta} \parallel_2^2 $$ where $\mathbf{\tilde{y}} = \mathbf{y} - \bar{y}\mathbf{1}_n$ Can someone tell me, what is the next step which will lead to a nice marginalization of $\mu$ and getting their result?


title: Correlation and Hotelling test

text: To find a Hotelling $T^2$ score it is necessary to calculate the covariance matrix and then invert it. Now, when the test is a two-sample $T^2$ test, the covariance matrix is a pooled matrix. When there is a high degree of correlation / collinearity in the variables the covariance could be close to a singularity (the determinant near zero). So it is recommended in statistics books to remove those highly correlated variables from the data matrix and to calculate the covariance matrix again until having a viable matrix that can be inverted. This is my question: **how do you determine the correlation in the variables when you have two groups?** My intuition tells me that you find the correlation matrix directly from your data matrix and you use this correlation matrix to determine which variables to eliminate. On the other hand, in a two-sample test you have a covariance matrix that is pooled: a weighted sum of two matrices. Should you find then the highly correlated variables from the covariance matrix of each one of the two groups and then eliminate them all from the data matrix? Which option is correct?


title: How to compare coefficients of a negative binomial regression for determining relative importance?

text: I'm working in R, using glm.nb (of the MASS package) to model count data with a negative binomial regression model. I'd like to compare the relative importance of each of my predictor variables regarding their impact on the response variable (note: the predictors each have quite different scales - sometimes by orders of magnitude). Unfortunately, the output from R gives me results as unstandardized ( _b_ ) coefficients ("estimates"). I'm hoping someone can give me a hint as to how to go about getting standardized ( _beta_ ) coefficients from the NB regression model... or another 'better' way to determine the relative importance of each of my predictors on my response variable. I've investigated several potential ways like:   1. using the R package 'relimpo' (as suggested in a comment to http://stats.stackexchange.com/a/7118), but it does not work on a NB regression model, thus completely changing the assumptions I should be accounting for and making the outcomes very different;    2. mean-centering and scaling my data, which changes the interpretation and makes it so that I can't use NB model due to response variables now having negative values;    3. scaling-only, so that I can still run a NB model... which I _thought_ would only affect the scale of the coefficients without changing their direction (viz., http://stats.stackexchange.com/a/29784 ) - but I do get some positive coefficients that flip to neg. and vice-verse... which seems strange to me and makes me wonder whether I'm making a mistake. I've benefited from looking at When should you center your data & when should you standardize? (and the suggested links from comments on the question such as http://andrewgelman.com/2009/07/when_to_standar/ and When and how to use standardized explanatory variables in linear regression and Variables are often adjusted (e.g. standardised) before making a model - when is this a good idea, and when is it a bad one?). Bottom line: I have not yet found a way to use a NB model in R (which I have statistically confirmed is more appropriate than lm, glm, or poisson for modeling my data) and still get at the relative importance - or at least to the standardized beta coefficients - for my predictors... The R scripts is something like this:               library("MASS")     nb = glm.nb(responseCountVar ~ predictor1 + predictor2 +        predictor3, data=myData, control=glm.control(maxit=125))     summary(nb)          scaled_nb = glm.nb(scale(responseCountVar, center = FALSE) ~ scale(predictor1, center = FALSE) + scale(predictor2, center = FALSE) +        scale(predictor3, center = FALSE), data=myData, control=glm.control(maxit=125))     summary(scaled_nb)


title: How to output training error when using cv.glmnet from the glmnet package in R?

text: I am currently using the `glmnet` package in R along with its' cross validation function `cv.glmnet`. As a reminder, `cv.glmnet` is an additional function that can be used to automatically run cross-validation on the regularization parameter, `lambda`. `cv.glmnet` can output a vector containing the mean cross-validation error for different values of `lambda`. That said, this is mean testing error, and I am wondering whether there is a way to have it output the mean training error as well.


title: Best practices for MCMC early stopping?

text: What are best and / or standard practices for MCMC early stopping? I have an algorithm which I want to compare with existing non-MCMC algorithms for accuracy and speed. When assessing the speed it's a bit tricky, since the speed is a function of the number of iterations I use in the Markov chain, and currently highly subjective. I'd like some kind of more objective way of deciding when to cut the chain; ideally some sort of 'best practice' that's applicable for reasonably well behaved Markov chains. Note that this is not some crazy algorithm that's going to jump between different apparently stable distributions every 2 weeks. It just goes up and down a bit, and then is stationary.


title: Regression coefficients without intercept

text: Could someone recommend a link or help me out here: where can I find the formula for the regression without an intercept, and how is it deriveed differently than the formula with the intercept? (matrix form)


title: Is there a reference that suggest using 30 as a large enough sample size?

text: I have read/heard many times that the sample size of at least 30 units is considered as "large sample" (normality assumptions of means usually approximately holds due to the CLT, ...). Therefore, in my experiments, I usually generate samples of 30 units. Can you please give me some reference which should be cited when using sample size 30?


title: Why is median age a better statistic than mean age?

text: * If you look at Wolfram Alpha ![alt text](http://i.stack.imgur.com/x0ILQ.png)   * Or this Wikipedia page List of countries by median age ![alt text](http://i.stack.imgur.com/QBtjT.png) Clearly **median** seems to be the statistic of choice when it comes to ages. I am not able to explain to myself why **arithmetic mean** would be a worse statistic. Why is it so? Originally posted here because I did not know this site existed.


title: Does the Bayesian posterior need to be a proper distribution?

text: I know that priors need not be proper and that the likelihood function does not integrate to 1 either. But does the posterior need to be a proper distribution? What are the implications if it is/is not?


title: Predicting class probabilities in regression based on area under the curve

text: Logistic regression models the log odds. That is for rv $Y$ which is binary logit$(Y=1)=X\beta$. Then with this model, you can estimate the class probabilities and hence prediction or classification is immediate. Ordinal regression does a similar thing with cumulative logits with which the class probabilities can then be calculated. Now looking at $AUC$ regression where Dodd and Peppe(2003) modeled the $AUC$, by $AUC$=$P(Y<Y^{\prime}$ $|$ $X,X^{\prime}$)=$f(X\beta)$ where $X$ are covariates. My question is how do we now do predictions with such a model.


title: Negative binomial distribution for bounded data

text: Why is the negbin distribution required when the analyzed count data is bounded? I don't really understand the following: > "The Poisson distribution can form the basis for some analyses of count data > and in this case Poisson regression may be used. This is a special case of > the class of generalized linear models which also contains specific forms of > model capable of using the binomial distribution (binomial regression, > logistic regression) or the negative binomial distribution where the > assumptions of the Poisson model are violated, in particular when the range > of count values is limited or when overdispersion is present." --- Wikipedia


title: Using multiple linear regression to distinguish two datasets

text: Consider two datasets, a study dataset with $n$ points and a control dataset with $n_c$ points, with $n$ I would like to evaluate the hypothesis that the study dataset has a different Y (in average or distribution) than that of the control dataset, after controlling for all independent variables $X_1$, $X_2$, $X_3$, $X_4$ simultaneously. Following a previous discussion, I applied multiple regression analysis to the two datasets. The coefficients of the linear regression are different, unsurprisingly. Since the control dataset is larger than the study one, I wanted to make sure that the difference was not the result of small(er) number statistics. So from the $n_c$ control observations I randomly selected a subset of $n$ and repeated the regression analysis, 10k times. The difference for one of the coefficients, the one with the largest value, is quite significant, at 2.7$\sigma$ when assuming a Gaussian distribution. Is this test conclusive in the sense that it proves that the datasets are different in what concerns Y? How would you suggest to do such a test? I played around with PCA but could not formulate the question in a concise fashion, but I am quite unhappy with the current dependence on the model assumption (linear).


title: R vs STATA which is more valuable to an actuary?

text: I'm going to school to be an actuary and I am just now starting to learn the basics of different statistical software packages. I have used both STATA and R, and I'm not sure which I like better. If there are any actuaries out there reading this what do you prefer? I currently use STATA for my econometrics class, and I have used R in my previous statistics class, and I must say I'm not sure which I like better. However I am very interested in improving my knowledge about whichever is most commonly used for the field I'm looking to go into. So any advice? And I appologize to anyone that doesn't like my question because of the topic but I don't know where else to post it to get real advice.


title: Computing composite score from variables that do not follow a likert scale

text: I am trying to develop a regression-weighted composite score. I have 4 variables a, b, c, d that are not necessarily linearly related that I would like to transform to A, B, C, D such that I can combine linearly by regression weights. So I am looking for a way of performing this transformation (a, b, c, d) --> (A, B, C, D). For ex: a, b, c, d may not be linearly related, but A = log2a,B = b, C = exp(c), D = 1/log(d) might be. How do I determine this relation? Can anyone guide me how to do this? I tried plotting scatter plots but they don't seem to convey much useful information. Should I just use the correlation? As a realtime example, say I want to combine Reputation score, #gold, #silver and #bronze badges and obtain a composite score.


title: can I adjust some confounding factors for a uncontrolled prospective observational study?

text: The patients with depression had a significant improvement in their severtity of disease after 8 weeks antidepressant treatment and I found the serum BDNF levels of the patients with depression after 8 weeks antidepressant treatment were lower compared to those in baseline. This study is a uncontrolled observational study, so we do not have a control group. Last month, a reviewer proposed a question that I should calculate the connection of the decrease of the BDNF and the course of treatment, and he thought using paired t test to analyze the data was not appropriate since there were many confounding factors needed to adjust such as length of disease, baseline severtity . I want to know if I can adjust the mentioned confouding factors when we do not have a control group and how can I calculate the connection of the BDNF decrease and the course of treatment?


title: Constrained kernel density estimation

text: Suppose you are trying to estimate the joint density $p(x,y)$ based on observed $(X,Y)$. However, you know that the marginal density $p(x)$ is uniform. How can you use this information to improve your density estimate?


title: PCA replicate data/index R

text: I have a set of date covering petrol prices. My example has two columns where each row represents a sequential date.                  unleaded diesel     1   1.39     1.35     2   1.3901   1.3502     3   1.3902   1.3501     .....      I generate eigen values: > my.eigen $values [1] 7.053791e-07 9.097811e-08               $vectors                  PC1        PC2     unleaded 0.6489256 -0.7608519     diesel   0.7608519  0.6489256      > my.eigen $values [1] 7.053791e-07 9.097811e-08 How can I produce a timeseries using the first, second or first two Components. That would give me something that replicates(but of course is different to) the original data? I would like to plot unleaded, diesel timeseries against the "replicated" data for different inclusion of components.


title: Comparison of before and after ordinal data across two groups

text: I would like to compare before-and-after "scores" on five quality of life questions across two groups. What I would like to know are:   1. Are the baseline scores across the two groups similar?   2. Is the change in scores significant (before and after) within each group?   3. Is the change in scores significant across both groups? Will an independent t-test suffice to answer question number 1? Will it be appropriate to perform a Wilcoxon signed rank test to answer question number 2? If not, what will be the more appropriate test? On question 3, is Mann-whitney appropriate? Thank you so much for the help.


title: In non-negative matrix factorization, does the first N eigenvector have N greatest variance?

text: I know for PCA, it's true that the first N eigenvectors have N greatest variance. But I'm not sure whether that's also true for NMF(Non-negative Matrix Factorization). For example, this method(Standard Nonnegative Matrix Factorization (NMF) [Lee2001], [Lee1999].): http://nimfa.biolab.si/nimfa.methods.factorization.nmf.html And are there articles that write about how to calculate the variance of first N eigenvectors? What is the best k(k means counts of eigenvector) for dimension reduction? Thanks!


title: Cross correlation for extremely highly correlated series

text: I have two time series showing very high dependency, in the order of 99.99% correlation, and I need to study their lead-lag relationship. So far I've been looking at the Pearson correlation with 7/8 decimals precision, but I wonder if there is a better approach. The problem is that most of the time the series value does not change: changes are concentrated in less than 1% of the observations. Should I extract a sub-sample based on the value changes in one of the series?


title: Kurtosis of a standardized Student's-t distribution?

text: I use the generalized form of the Student's-t distribution: \begin{align*} f(l|\nu ,\mu ,\beta) = \frac{\Gamma (\frac{\nu+1}{2})}{\Gamma (\frac{\nu}{2}) \sqrt{\pi \nu} \beta} \left(1+\frac{1}{\nu}\left(\frac{l - \mu}{\beta}\right)^2 \right)^{\text{$-\frac{1+\nu}{2}$}} \end{align*} I want to have standardized version, i.e. mean zero and a variance of one. Therefore, I set $\mu=0$ and \begin{align*} \beta=\sqrt{\frac{\nu-2}{\nu}} \end{align*} which ensures, that the variance is equal to one. If I now insert this I get after some derivations \begin{align*} f(l|\nu) =(\pi (\nu-2))^{-\frac{1}{2}}\Gamma \left(\frac{\nu}{2} \right)^{-1} \Gamma \left(\frac{\nu+1}{2} \right) \left(1+\frac{l^2}{\nu-2} \right)^{-\frac{1+\nu}{2}} \end{align*} Now my question is: What is the formula for the kurtosis? Is it still $\frac{6}{\nu-4}$? E.g. consider these data and the following R code:               pinumber<-3.141592653589793     startvalue<-2          loglikstandardizedt <-function(par){     if(par>0) return(-sum(log((pinumber*(par-2))^(-1/2)*gamma(par/2)^(-1)*gamma((par+1)/2)*(1+standresidsapewma^2/(par-2))^(-(1+par)/2))))     else return(Inf)     }          optim(startvalue, fn=loglikstandardizedt, method="Brent",lower=2,upper=250)     param = optim(startvalue,loglikstandardizedt, method="BFGS")$par      If I look at the plot, to see how good the fit is, I do the following code:                # control output         denstiystandtresid<-function (x) (pinumber*(param-2))^(-1/2)*gamma(param/2)^(-1)*gamma((param+1)/2)*(1+x^2/(param-2))^(-(1+param)/2)              plot(density(standresidsapewma),ylim=c(0,0.8))         curve(denstiystandtresid,col="red",add=TRUE)      This gives me the following plot: ![plotgraph](http://i.stack.imgur.com/bRJiE.png) As you can see, the fit is, let's say fairly ok. Now, I am interested in the excess kurtosis. The data has the excess kurtosis of               kurtosis(standresidsapewma)      which gives `0.6470055` I would expect, since the fit is quite ok in the tails, that the fitted distribution has almost the same excess kurtosis, but if I calculated it via the following way (the estimate output for $\nu$ is 8.85009): $\frac{6}{\nu-4}=\frac{6}{8.85009-4}=1.23709$? which is pretty much more than 0.64. This seems to be wrong to me, since I believe, that the fit in the tail is quite ok, so the kurtosis should be almost the same? Is my formula for calculating the ex kurtosis in case of a standardized Student's-t distribution wrong? Or what is my mistake?


title: Fisher's exact test vs chi-squared with infection data

text: I have data from an infection study that I did and I'm trying to figure out if the differences I have observed are significant or not. I have infections with 4 different bacterial strains, with 10 mice infected per strain (40 mice total). Then a colleague blindly scores the inflammation induced by the bacteria on a scale of inflamed, mildly inflamed, and not inflamed. I understand that this creates data that is qualitative and ordinal, since the inflammation is ranked. I'm trying to figure out if I should be using a Fisher's exact test or a chi- squared test (or something else?), as well as what the minimum group size should be. I already observe big differences comparing the wild-type infected mice to my mutant infected mice with only 10 mice per group. When I run a chi-squared test in Graphpad it tells me that there are significant differences but the calculations are only valid when all expected values are greater than 1.0 and at least 20% of the expected values are greater than 5. Can anyone suggest what to do? I'd like to use the least number of mice possible.


title: Problem with Spearman correlation in case of many ties

text: I'm using the Spearman correlation coefficient test for multiple non- parametric ranked categories in excel. I used the built in rank function in excel to calculate the rank of each individuals but i have many ties because there are 10 to 208 categories for 2k instances. In the variable with 10 categories excel ranks go directly from 1 to 2100 with no intermediary ranks. Should it be doing this? The spearmens test is giving high r values like .25 for 2747 instances and very low p values for the two ranks .25 but pearsons on the two ranks gives negligible r values like .05. This is the formula I followed: http://blog.excelmasterseries.com/2010/09/spearman-correlation-coefficient- test.html Thank you


title: What is the unique thing about super-population for a meta-analysis?

text: The Encyclopedia of survey research indicates that super-population is similar to what may be termed as infinite population. However, in a meta-analysis of effect-sizes, we use outcomes of sample studies that seem to originate from the specific target populations i.e, we may treat the study outcomes from a single specified population (e.g. studies on US population). It is not clear in terms of sampling and meta-analysis, what is the key assumption in case of random-effects model(for a meta-analysis) proposed by Hedges and Olkin (1985).


title: How to assess stability of daily time series in sentiment analysis?

text: I developed a measure of "sentiment" and I have time based data and used the measure to derive a daily sentiment time series. I am looking for some way to establish reliability or maybe stability. For example, when there are dramatic changes in the "sentiment" is this "real". Two ideas for doing this numerically start with re-sampling the data. By this I mean, drawing random samples ( on a per day basis so that for example if I have 100 data points per day, I would have 25 data points per day in each of the random samples). Then, one way of establishing reliability would be testing the average variance in the resulting curves (distributions). The other way was to take each of the resulting data sets and perform pairwise the Kolmogorov-Smirnov test was and check the significance values. Can anyone comment on either of these approaches or suggest another?


title: In a matched case-control study, what do I use for hypothesis testing in the descriptive statistics or the "univariate" associations?

text: I am doing a case-control study with 80 disease cases matched 1:3 to non- diseased controls and examining whether they had a binary exposure prior to developing the disease. I am using multivariable conditional logistic regression for the main analysis. My question is, for the "descriptive" or "univariate" table where we try to investigate simple associations of each (non-matched) variable with the outcome, do I just use something like a chi- square test and ignore the matching, or do I use main effects conditional logistic regression (or something else)? Any help would be greatly appreciated.


title: How can I transform Ornstein-Uhlenbeck model into state space form?

text: As I say in the subject, How can I put the model $d x_t = \eta\, (\overline{x} - x_t)\,d t + \sigma\, x_t\,d W_t$ into state space form? I mean, which are the observation and transition matrices?


title: How do I calculate the variance of the OLS estimator $\beta_0$, conditional on $x_1, \ldots , x_n$?

text: I know that $$\hat{\beta_0}=\bar{y}-\hat{\beta_1}\bar{x}$$ and this is how far I got when I calculated the variance: \begin{align*} Var(\hat{\beta_0}) &= Var(\bar{y} - \hat{\beta_1}\bar{x}) \\ &= Var((-\bar{x})\hat{\beta_1}+\bar{y}) \\ &= Var((-\bar{x})\hat{\beta_1})+Var(\bar{y}) \\ &= (-\bar{x})^2 Var(\hat{\beta_1}) + 0 \\ &= (\bar{x})^2 Var(\hat{\beta_1}) + 0 \\ &= \frac{\sigma^2 (\bar{x})^2}{\displaystyle\sum\limits_{i=1}^n (x_i - \bar{x})^2} \end{align*} but that's far as I got. The final formula I'm trying to calculate is \begin{align*} Var(\hat{\beta_0}) &= \frac{\sigma^2 n^{-1}\displaystyle\sum\limits_{i=1}^n x_i^2}{\displaystyle\sum\limits_{i=1}^n (x_i - \bar{x})^2} \end{align*} I'm not sure how to get $$(\bar{x})^2 = \frac{1}{n}\displaystyle\sum\limits_{i=1}^n x_i^2$$ assuming my math is correct up to there. Is this the right path? \begin{align} (\bar{x})^2 &= \left(\frac{1}{n}\displaystyle\sum\limits_{i=1}^n x_i\right)^2 \\ &= \frac{1}{n^2} \left(\displaystyle\sum\limits_{i=1}^n x_i\right)^2 \end{align} I'm sure it's simple, so the answer can wait for a bit if someone has a hint to push me in the right direction.


title: Fitting distributions in R

text: I hope you are very well. I have a big dataset (~9 million registries) and I have 2 variables $X$=purchase amount and $Y$=frequency of purchase. I would like to know what distribution should I use for each variable and fitting its distribution with `R`. For example, for $X$ I'm testing with package 'mixtools' but the results are not really good. I want to create a simulation model, that simulates that variables. First of all, I'm selecting appropriate probability distributions (and associated parameters) to describe the behavior of each uncertain input variable. Thanks. UPDATE: I have replaced the pdf with a histogram for each variable. ![](http://i.imgur.com/WTcGoQZ.jpg?1)


title: How can I change my fitting method to force better relative accuracy instead of just absolute accuracy?

text: I have some data which I want to use for interpolation and extrapolation. I would (very much) prefer a function that I can plug values in instead of using cubic splines or some such technique so I resorted to least squares fitting. Here are some plots with explanation interspersed. ![enter image description here](http://i.stack.imgur.com/mhgie.jpg) In this first picture above, I have four of my good data segments. Here are my goals. I want to use an appropriate form for the fitting function and then use (most likely) non-linear least squares fit to determine the parameters because I suspect the fitting function will be non-linear in its parameters-to-be- estimated. I want to place "heavy emphasis" if necessary, on the first half of the segment meaning between 0.8 and 7-8mHz-ish. I want the fit to be much more accurate for smaller frequencies. Don't care much for the higher frequencies. This leads into my second goal because I then want to use my fit for extrapolation between 0mHz and 0.8mHz. This is just to show the overall picture. So seeing the bump, my first attempt was to try a gaussian plus a linear term $$y = Ae^{-(\frac{x-B}{C})^2}+Dx+E$$ the problem with that was that the slope of the line is too large and the values-to-be-extrapolated were negative for example at $x$=0.5mHz. In addition, for some of the data like the red curve, the interpolated values like at 1mHz were negative. And this is a big no-no because all values must be positive. Then my next attempt is a gaussian plus power law $$y = Ae^{-(\frac{x-B}{C})^2}+Dx^E$$ so all of the fitted curves from now on, including the one in this first picture are of the form gaussian plus power law. This first picture includes food data segments with their fits as well as the correlation coefficients. ![enter image description here](http://i.stack.imgur.com/k7lV1.jpg) This second picture is just the first one zoomed in to show that even though the fits are "good", the are rather "bad" in this range. Least squares goodness is measured in terms of absolute difference which in this case is small. But I am also concerned with relative difference. In this case for the first point on the red curve, you can see that my interpolated value is four times larger then the actual measured value. ![enter image description here](http://i.stack.imgur.com/Mr9xh.jpg) Here I show the log-log plot of the data with their fits just to emphasize the difference. The fits are "way-off" in log-log space and I need just the opposite, something like the fits going through the first three points or something. There is a clear linear trend in this log-log plot, hence the power law I added but I was hoping that this linear trend would just nicely continue all the way down to zero and I can use it for extrapolation. I also see the intersection of the top two fits as a "problem". Straight lines for the top two curves just decreasing to the left shouldn't cross. Then my next attempt was to use the same form, gaussian plus a power law, but to truncate half of my data from the right and repeat all this to get a better fit. ![enter image description here](http://i.stack.imgur.com/VMvcI.jpg) Here I only consider the first half of the data, from 0.8mHz to 8.1mHz and then reran the fits. The above picture shows the results. ![enter image description here](http://i.stack.imgur.com/QNkQ7.jpg) Again, just to show what's happening, I zoom in. The fits are "better" absolutely and relatively but the relative errors are still too "large". Now the interpolated value is only twice as large as the measured value. ![enter image description here](http://i.stack.imgur.com/rVNaJ.jpg) This above picture shows the same truncated data and fits on those truncated data on a log-log scale and we still see that crossing in the red ellipse. So now here are my questions. And I know these are probably "too many questions for one post" but they are all related and the answer to one may alter the answer to another.   1. Looking at all these plots and log-log plots, is there another more appropriate fitting form? I saw a straight-line and a bump in log-log space so used gaussian plus power law.   2. All of these fits are being done to the original data without any kind of a transformation. I would love to just take log of the data and then fit it and it would work wonderfully (because I tried that) but that is of course wrong with least squares. Is there something I can do with least squares where relative errors will be taken into account instead of just absolute errors? Like mess around with my fitting form? Like replace $x$ with $\log(x)$ and $y$ with $\log(y)$ and then fit with least squares?   3. Is there another technique, analogous to least squares, besides least squares, which I should look at? Something non-parametric? I don't know much about stats/regression in general but something non-parametric returning a nice function form I can use for the entire interval (truncated or not) would be ideal (instead of 10-piecewise functions like how you get from a cubic-splines). Some technique which minimizes sum of the squares in a different way or even minimizes something else entirely in a non-parametric way. Something where I provide the data and the fitting form, etc., and it returns parameter-estimates. Because non-parametric, I can use in log-log space, right? I also tried weighted least squares but it just gets ridiculous. I need to weigh the first three points like 30,000 times more then the last 17 points and it still doesn't quite work. The relative error is still too large for smaller frequencies and it seriously screws up the fit for the larger frequencies (for obvious reasons) so the fit sucks for large AND small frequencies. The following is an example of a picture of what I am trying to achieve here, with a different data set, on which gaussian+linear worked. ![enter image description here](http://i.stack.imgur.com/g5L6y.jpg) Top left point, first in red, the relative error (fit/measured) is 1.07. This example had much better fits because the absolute errors were larger for smaller frequencies so minimizing them was a priority. With my data I am currently working with, the trend is reversed so smaller frequencies are "ignored" and are "not-as-important" so least squares focuses on higher frequencies which I don't really care about anyway. Help! Thank you, I know this got long but wanted to explain what I am trying to achieve and what I have tried.


title: Statistics of sample correlation for uniformly distributed samples

text: I am computing the sample correlation between two vectors of uncorrelated and uniformly distributed samples using MATLAB. More precisely, I compute $$ r_N=\frac{1}{N}\sum_{i=1}^N x_{i}\, y_{i}, $$ where $E(X)=E(Y)={\rm cov}(X,Y)=0$, and $x_{i}, y_{i}$ are drawn from a uniform distribution in $[-A/2,A/2]$. It is clear that $\lim_{N \rightarrow \infty} r_N=0$. However, I would like to compute an estimate of $\tilde{N}$ such that, given $\epsilon >0$ and $\alpha >0$, we have $ {\rm Pr}(|r_N| < \epsilon) < \alpha$ for all $ N > \tilde{N}$. In other words, I need an expression for the sample statistic $r_N$. I formulated this question for uniformly distributed data; however I'm interested to know what happens for other distributions (say normal). Any ideas?


title: Fitting a gaussian mixture model using stochastic gradient descent

text: I'm working on an online category learning model which uses stochastic gradient descent to fit a gaussian mixture model. The model is based on the online learning model used in Toscano & McMurray (2010). While gradient descent seems to be working fairly well to estimate the means and frequencies/mixing probabilities of the categories, I'm having issues with estimating the covariances of mixture components. The partial derivatives I've been using for the gradient descent update come from Petersen & Pedersen (2008) (p. 44) Starting with $p(x) = \sum _k \rho_k \mathcal N_x(\mu_k,\Sigma_k)$ Petersen & Pedersen give the partial derivative with respect to the covariance matrix $\Sigma$ as $\frac{\delta \ln p(x)}{\delta \Sigma_j}=\frac{\rho_j\mathcal N_x(\mu_j,\Sigma_j)}{\sum _k \rho_k \mathcal N_x(\mu_k,\Sigma_k)}\frac{1}{2}[-\Sigma_j^{-1}+\Sigma_j^{-1}(x-\mu_j)(x-\mu_j)^T\Sigma_j^{-1}]$ The gradient descent step for each $\Sigma_j$, as I've got it implemented in Python is (this is a slight simplification and the $\Delta\Sigma$ for all components is calculated before performing the update):               j.sigma += learning_rate*(G(x)/M(x))*0.5*(-inv(j.sigma) + inv(j.sigma).dot((x-j.mu).dot((x-j.mu).transpose())).dot(inv(j.sigma)))      Where j is an object representing the $j$th component of the mixture and j.sigma and j.mu are that component's mean and variance. G(x)/M(x) is standing in for some code that calculates $\frac{\rho_j\mathcal N_x(\mu_j,\Sigma_j)}{\sum _k \rho_k \mathcal N_x(\mu_k,\Sigma_k)}$ So, I'm wondering if there's something wrong with my code (highly likely) or if this is just a really bad way to fit this kind of model when dealing with data with more than two dimensions (See Toscano & McMurray for algorithms for univariate and bivariate data that definitely work). references: Toscano, J. C., & McMurray, B. (2010). Cue integration with categories: Weighting acoustic cues in speech using unsupervised learning and distributional statistics. Cognitive Science, 34, 434-464. Petersen & Pederson. The Matrix Cookbook, Version: November 14, 2008


title: Calculating threshold turn value using mean centered quadratic terms in regression

text: I have a question on calculating the threshold value or value at which the quadratic relationship turns. The formula for calculating the turn is at $x = -b/2a$, following from $ax^{2}+bx+c$. My question is this: when using the mean centered quadratic terms, do you add the mean value back to calculate the threshold turn value on the non-centered term (for purposes of interpretation when writing up results and findings).


title: Statistical significance of genes associated with multiple phenotypes

text: I have a dataset of gene-phenotype association in this format. I am looking at some combination of phenotypes and genes shared between combination. I would like to use a statistical test to show that the genes shared between two phenotypes are statistically significant using a p-value or a similar measure. For example:               22 genes are associated with Phenotype1      205 genes are associated with Phenotype2      9 genes are common between two phenotypes         I want to assess whether the number of genes common to two phenotypes are statistically significant or just a random observation. _I have phenotype information for 4035 genes; I assume that human genome contains 42, 071 genes_ How do you address this problem (preferably in R), what statistical test you would recommend and why ? PS. I have asked this question previously at BioStar


I am sorry for a more basic question, but I was unable to find any good sources on google. How is a bifactor IRT model different from a factor analysis? How would you describe their key differences? Any references you could direct me to would be very helpful. Thank you for sharing your knowledge!


title: Distribution family of the mean of iid random variables

text: Let $X_1,X_2,...,X_n$ be $n$ iid random variables following some distribution family $A$ (beta, gamma, normal, etc.). Does the mean $\bar{X}=\frac{1}{n}\sum{X_i}$ also follow a $A$-like distribution (with different parameters)? I know it is approximately normal due to the Central Limit Theorem, but does it follow the same family $A$? I'm particularly interested in Beta distributions. Thanks


title: Unable to fit negative binomial regression in R (attempting to replicate published results)

text: Attempting to replicate the results from the recently published article, > Aghion, Philippe, John Van Reenen, and Luigi Zingales. 2013. "Innovation and > Institutional Ownership." American Economic Review, 103(1): 277-304. (Data and stata code is available at http://www.aeaweb.org/aer/data/feb2013/20100973_data.zip). Am having no problem recreating the first 5 regressions in R (using OLS and poisson methods), but am simply unable to recreate their negative binomial regression results in R, while in stata the regression works fine. Specifically, here's the R code I've written, which fails to run a negative binomial regression on the data:               library(foreign)     library(MASS)     data.AVRZ <- read.dta("results_data2011.dta",                       convert.underscore=TRUE)     sicDummies <- grep("Isic4", names(data.AVRZ), value=TRUE)     yearDummies <- grep("Iyear", names(data.AVRZ), value=TRUE)     data.column.6 <- subset(data.AVRZ, select = c("cites",                                      "instit.percown",                                      "lk.l",                                      "lsal",                                      sicDummies,                                      yearDummies))     data.column.6 <- na.omit(data.column.6)          glm.nb(cites ~ .,        data = data.column.6,        link = log,        control=glm.control(trace=10,maxit=100))      Running the above in R, I get the following output:               Initial fit:     Deviance = 1137144 Iterations - 1      Deviance = 775272.3 Iterations - 2      Deviance = 725150.7 Iterations - 3      Deviance = 722911.3 Iterations - 4      Deviance = 722883.9 Iterations - 5      Deviance = 722883.3 Iterations - 6      Deviance = 722883.3 Iterations - 7      theta.ml: iter 0 'theta = 0.000040'     theta.ml: iter1 theta =7.99248e-05     Initial value for 'theta': 0.000080     Deviance = 24931694 Iterations - 1      Deviance = NaN Iterations - 2      Step halved: new deviance = 491946.5      Error in glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,  :      NA/NaN/Inf in 'x'     In addition: Warning message:     step size truncated due to divergence      Have tried using a number of different initial values for theta, as well as varying the maximum number of iterations with no luck. The authors' supplied stata code works fine, but I still can't seem to coerce R into making the model work. Are there alternative fitting methods for glm.nb() that may be more robust to the problem I'm encountering?


title: Shall I use sqrt of J-S divergence or K-L divergence as a measure of deviation from uniform pdf?

text: Some colleagues are using 'normalized' K-L divergence to measure deviation of unit-area histogram (pseudo discrete pdf) from corresponding uniform distribution ($N$ equal-length bins). Maybe it is better to use true metric like sqrt of J-S?