title: How to combine the forecasts when the response variable in forecasting models was different?

text: **Introduction** In forecasts combination one of the popular solutions is based on the application of some information criterion. Taking for example Akaike criterion $AIC_j$ estimated for the model $j$, one could compute the differences of $AIC_j$ from $AIC^* = \min_j{AIC_j}$ and then $RP_j = e^{(AIC^*-AIC_j)/2}$ could be interpreted as the relative probability of model $j$ to be the true one. The weights then are defined as $$w_j = \frac{RP_j}{\sum_j RP_j}$$ **Problem** A difficulty that I try to overcome is that the models are estimated on the differently transformed response (endogenous) variables. For example, some models are based on annual growth rates, another - on quarter-to-quarter growth rates. Thus the extracted $AIC_j$ values are not directly comparable. **Tried solution** Since all that matters is the difference of $AIC$s one could take the base model's $AIC$ (for example I tried to extract `lm(y~-1)` the model without any parameters) that is invariant to the response variable transformations and then compare the differences between the $j$th model and the base model $AIC$. Here however it seems the weak point remains - the difference _is_ affected by the transformation of the response variable. **Concluding remarks** Note, the option like "estimate all the models on the same response variables" is possible, but very time consuming. I would like to search for the quick "cure" before going to the painful decision if there is no other way to resolve the problem.


title: R: multiple linear regression model and prediction model

text: Starting from a linear `model1 = lm(temp~alt+sdist)` i need to develop a prediction model, where new data will come in hand and predictions about `temp` will be made. I have tried doing something like this: `model2 = predict.lm(model1, newdata=newdataset)` However, i am not sure this is the right way. What i would like to know here is, if this is the right way to go in order to make prediction about `temp`. Also i am a bit confused when it comes to the `newdataset`. Which values should be filled in etc.? Any help is welcome.


What model to use to predict behavior of a group based on past behavior of group members. I am using an analogy of sports team. Say we have data for individual players of sports. The data is the results of games between 2 teams and the members of each team. From this we can see how many times players have played together and the results each time. Win or Lose. Today we know the roster of 2 teams playing. Is there a model to use that will take into account the past performance of each player when they have played with the other players on the field today and use it to predict the winner of the game?


title: Binary features for prediction

text: I have a set of relatively long ($\sim 1000$) binary features with scalar values $[0-10]$ attached to them. My aim is to write a predictor that learns to map the features to the $[0-10]$ interval to predict new features when given a new binary vector. I used SVM and Lasso with leave-one-out performance analysis, but both always end up predicting the mean value of the distribution (correlates to the histogram of all the feature - scalar distribution). The histograms are also rather norm / Rayleigh distributions. Suggestions for algorithms / feature space mapping? My main problem is that I am dealing with binary features for the first time.


title: How do economists quantify black market operations?

text: I was doing a lot of research into organized crime in East Asia for a project as a favor to an author-friend of mine, and I noticed that there were noted economists and journalists who, in conjunction, estimate the value of black market operations worldwide. What is the methodology for this? How do these numbers come together? Are these numbers captured in traditional market models?


title: Boruta test and naive bayes classification

text: I am currently using Boruta to test which feature is the most important to be used in my model development. For example, I have 3 features(X,Y,Z).Boruta test give the highest importance is Z. However when I test the feature with naive Bayes classifier, feature X is better performance. May I know why and how this situation happen? Thanks a lot


title: Assumption of Gaussian distribution of acceleration

text: I have a data set consisting of noisy position values of a trajectory of a human hand. I want to estimate a generative model of these trajectories, and the obvious choice is a Kalman Filter/linear dynamical system. Since the relationship of position/velocity/acceleration is well known, I chose a LDS with 3 latent variables (for each dimension) and hard coded the transition and emission matrices. I tuned the variances as diagonal covariances via "occular inspection". (I maybe should have learned all the parameters via EM, with possibly more stats, with a prior on the model parameters somewhat, I know.) Having done that, I realized that the LDS struggled with sudden changes of direction. Thus I wondered whether the assumption that acceleration is actually Gaussian distributed is correct. Now I have several questions.   1. How would I check whether the Gaussian assumption is a good one? One problem is, that I cannot get a reliable sample from the accelerations, since even the position signals have outliers, and so the outliers of the second derivative will be even more severe.)   2. What state space models are there that do not assume Gaussianity all over the place? Which of those have inference fast enough to be usable in a real time system with ~20 predictions a second?   3. Does anyone have any experiences with similar settings and might have any advice on which models to use?


title: Question about cointegration

text: I was reading about cointegration (page 5) and it said > "If a linear combination of combination of a set of $I(1)$ variables is > $I(0)$, then the variables are cointegrated." Does this mean that any 2 $I(1)$ variables, no matter what they are will definitely be cointegrated? If yes, can you give the proof. If no, could you give a counterexample. I'm very curious about this topic.


title: What does $\beta$ tell us in linear regression analysis?

text: I came across this note from a book: "..the correlation between the active and passive portfolios is greater when the $\beta$ of the active portfolio is higher.." The author runs regression of active portfolio on passive portfolio. This does not make sense to me. If I understand correctly, beta in a single variable linear regression analysis gives us the best relationship/fit between two variables (or the slope of independent vector on dependent vector), but it is not indicative of strength of that relationship. We need to look at correlation or $R^2$ value to determine how good is the relationship between the variables. As long as beta is statistically significant, we can look at $R^2$ to determine the strength of relationship. Here are couple of examples:   1. Consider three variables $A_1$, $A_2$ (dependent variables) and $B$ (independent variable). If $A_1$ and $B$ have $\beta=2.0$; correlation=0.8, and $A_2$ and $B$ have $\beta=0.5$; correlation=0.8, then both $A_1$ and $A_2$ are equally well explained by $B$.   2. Consider two series: A=$\\{1,2,3,4,5,6,7,8,9,10\\}$ and $B=2\times A$. In this case, correlation between A and B is 1. However, regression of $A$ on $B$ gives $\beta$ as 0.5 and $B$ on $A$ gives 2.0. However, each variable is completely explained by the other. Is my understanding correct? Please highlight if I am missing something here.


title: Non-linear mixed effects regression in R

text: Surprisingly, I was unable to find an answer to the following question using Google: I have some biological data from several individuals that show a roughly sigmoid growth behaviour in time. Thus, I wish to model it using a standard logistic growth               P(t) = k*p0*exp(r*t) / (k+p0*(exp(r*t)-1))      with p0 being the starting value at t=0, k being the asymptotic limit at t->infinity and r being the growth speed. As far as I can see, I can easily model this using nls (lack of understanding on my part: why can I not model something similar using standard logit regression by scaling time and data? EDIT: Thanks Nick, apparently people do it e.g. for proportions, but rarely http://www.stata-journal.com/article.html?article=st0147 . Next question on this tangent would be if the model can possibly handle outliers >1). Now I wish to allow some fixed (mainly categorical) and some random (an individual ID and possibly also a study ID) effects on the three parameters k, p0 and r. Is nlme the best way of doing this? The SSlogis model seems sensible for what I am trying to do, is that correct? Is either of the following a sensible model to begin with? I cannot seem to get the starting values right and update() only seems to work for random effects, not fixed ones - any hints?               nlme(y ~ k*p0*exp(r*t) / (k+p0*(exp(r*t)-1)), ## not working at all (bad numerical properties?)                 data = data,                 fixed = k + p0 + r ~ var1 + var2,                 random = k + p0 + r ~ 1|UID,                 start = c(p0=1, k=100, r=1))          nlme(y ~ SSlogis(t, Asym, xmid, scal), ## not working, as start= is inappropriate                 data = data,                 fixed = Asym + xmid + scal ~ var1 + var2, ## works fine with ~ 1                 random = Asym + xmid + scal ~ 1|UID,                 start = getInitial(y ~ SSlogis(Dauer, Asym, xmid, scal), data = data))      As I am new to non-linear mixed models in particular and non-linear models in general, I would appreciate some reading recommendations or links to tutorials / FAQs with newbie questions.


title: Slope of line when Correlation = 0 is not always 0?

text: Practicing for a statistics final in 2 days. When is the following true?               The correlation between the variables x and y is zero.          However, the best curve that relates x to y is not a horizontal line.      Any ideas? Thank you.


title: Un-transformable Seemingly Continuous Data? If not, What to Do With Non-Parametric Test Results - What is the Next Action?

text: Note: If you want to skip to the summary of my 3 main questions, they are at the very bottom. I am attempting to produce a binary logistic regressional model for determining whether or not the data I am given is correct. I have a particular column of data that is giving me quite a bit of trouble. The column of data is very skewed. It looks like a very, very sharp exponential decay. It is used in conjunction with other indepedent variables that are linear. In this case I am wanting to know when a zip code I am given from a certain data source with X categorical input is correct with at least 90% confidence (ideally higher). I already know ~25% of the zip codes are correct - as I've compared them to data I know to be true. Keep in mind this is before scrubbing as I'll explain shortly. After scrubbing ~40% of zip codes are true. In the future when I have no true data to compare the returned dataset with, I'll know how many zip codes are correct, but I won't know which ones to pick without figuring this (among a number of others problems) out. I am given a lot of other data within the dataset I am returned as well that I've considered trying to somehow proportionate (or do something) with to correct for this problem, but I feel like doing so is the wrong action. This would then cause dependence between my supposedly independent variables - and would thus be pointless. This is correct?...right? Actually, barring a counterexample that generally seems obvious to me. I am somewhat new to statistics, but I am learning as quickly as possible so forgive me for any false assumptions, etc. I make here. I have software that has allowed me to apply nearly every non-normalized data transformation technique you can think of to this - but the AD p-value always comes back with a value much less than .05. About 10+ transformations were attempted. I list them below the dataset. I am pretty sure I am going to have to go the non-parametric route, which I generally understand - except for one key part: Everything I can find about non-parametric is about "tests." Like the Mann Whitley test, etc. I can't find anywhere or way I can actually transform the data non-parametrically. **I can't figure out what to do with the result of my non-parametric "tests". What do you do with them? What is the action taken from the results of the test?** The only thing I can see that I could do here is straight-up using the data from this table, figuring out a proportion based upon how likely it is correlated with being correct for each number, assuming that number is categorical rather than ordinal at first (and definitely not continuous), and then sorting that list of proportions in a ordered fashion according to the likelyhood they dictate that the zip code given is the correct one. I'd come up with a list of ordinal variables I could use. The huge problem is this doesn't account for standard deviations, number of times the original variable was listed, ect.. That is bad, so pending further enhancement of suggestions that strategy is out the window. **Another problem: Blanks** I know that you are not supposed to delete blanks as blanks mean something if they are there - but even after inquiring the source of the returned data I get, the only answers I receive are simply "we just didn't have the information." I had no idea what to do there and considered trying to use regressional models in an attempt to fill in the blanks. However, since 95% of the rows where a blank appeared in one of the columns equated to the given zip code being wrong ~98% of the time anyway...I did the unthinkable: I just decided to delete them and add this "scrubbing" into my formula, which would essentially make it a simple algorithm. I assumed as long as did the same scrubbing technique everytime, the eventual formulaic model I came up with would work the same as future data would equivalently be scrubbed of blanks. This is assuming the way the source of data is derived remains the same every time. If anyone has a better idea, suggestion, or correction for this, let me know. Here is the data. I hope the list isn't too long - and yes it always begins at 1 and the highest value it ever reaches is 999. It's beautifully skewed at both ends. This has been scrubbed of rows of data that contain blanks. The blanks do not only appear in the Accuracy Radius column though - this happens only about .02% (so yes .0002) of the time. The blanks could appear in a variety of columns as well - all indicating the zip code returned is incorrect ~98% of the time as previously stated. I realize this may potentially affect the result I get when attempting to transform this specific column of data in some way, but if I perform the said scrubbing every time as I mentioned I theorize it's fine. If anyone would like the data with blanks, let me know. Here is the data and below it are my summarized questions. Obviously column four would just be the sum of the previous two columns - it's the number of times that specific accuracy radius number was listed. I just listed it for convenience. I hope I formatted it correctly - it was quite difficult to do.               Accuracy Radius (km)    #Zip=Correct    #Zip=Incorrect      #Rows AccR is Listed                                                                     1               1041                    327                 1368             2               764                     727                 1491             3               681                     1070                1751             4               527                     760                 1287             5               339                     605                 944             6               248                     405                 653             7               157                     320                 477             8               108                     245                 353             9               87                      174                 261             10              70                      124                 194             11              78                      105                 183             12              47                      103                 150             13              41                      73                  114             14              46                      56                  102             15              37                      59                  96             16              23                      50                  73             17              31                      61                  92             18              15                      47                  62             19              27                      36                  63             20              25                      41                  66             21              19                      30                  49             22              15                      43                  58             23              18                      52                  70             24              9                       47                  56             25              10                      24                  34             26              11                      30                  41             27              7                       32                  39             28              12                      25                  37             29              16                      32                  48             30              10                      31                  41             31              9                       23                  32             32              4                       25                  29             33              11                      24                  35             34              10                      18                  28             35              3                       21                  24             36              7                       20                  27             37              12                      17                  29             38              6                       9                   15             39              5                       13                  18             40              5                       18                  23             41              2                       19                  21             42              5                       24                  29             43              2                       16                  18             44              3                       12                  15             45              7                       12                  19             46              2                       21                  23             47              1                       9                   10             48              2                       15                  17             49              7                       12                  19             50              5                       11                  16             51              7                       13                  20             52              3                       15                  18             53              1                       12                  13             54              4                       6                   10             55              3                       12                  15             56              4                       13                  17             57              4                       5                   9             58              5                       12                  17             59              4                       10                  14             60              4                       6                   10             62              2                       7                   9             63              1                       8                   9             64              3                       9                   12             65              0                       6                   6             66              2                       5                   7             67              2                       13                  15             68              2                       14                  16             69              1                       4                   5             70              2                       7                   9             71              2                       12                  14             72              1                       8                   9             73              0                       9                   9             74              1                       9                   10             75              3                       2                   5             76              1                       8                   9             77              1                       8                   9             79              5                       8                   13             80              2                       5                   7             81              5                       4                   9             82              0                       12                  12             83              0                       5                   5             84              3                       5                   8             85              0                       5                   5             86              2                       3                   5             87              9                       8                   17             88              0                       8                   8             89              1                       8                   9             90              2                       8                   10             91              0                       10                  10             92              0                       4                   4             93              0                       10                  10             94              2                       5                   7             95              2                       6                   8             96              1                       5                   6             97              2                       7                   9             98              0                       9                   9             99              1                       7                   8             100             2                       11                  13             102             0                       3                   3             103             2                       5                   7             104             1                       7                   8             105             0                       10                  10             106             0                       5                   5             108             1                       6                   7             110             2                       4                   6             111             1                       9                   10             114             2                       3                   5             115             0                       4                   4             116             0                       5                   5             117             2                       3                   5             119             1                       1                   2             120             1                       3                   4             122             2                       5                   7             125             1                       4                   5             126             0                       2                   2             127             0                       4                   4             128             0                       1                   1             129             1                       2                   3             130             0                       3                   3             131             0                       3                   3             132             1                       3                   4             133             0                       6                   6             134             0                       5                   5             137             0                       3                   3             139             0                       3                   3             141             2                       8                   10             143             0                       3                   3             145             2                       2                   4             146             0                       4                   4             149             1                       1                   2             150             0                       7                   7             152             0                       6                   6             153             0                       6                   6             156             1                       4                   5             157             1                       3                   4             160             1                       3                   4             161             0                       17                  17             163             0                       1                   1             165             0                       1                   1             166             0                       5                   5             170             0                       6                   6             171             0                       4                   4             173             1                       13                  14             178             0                       4                   4             180             0                       18                  18             184             0                       15                  15             185             0                       6                   6             186             0                       1                   1             187             0                       2                   2             193             0                       2                   2             197             0                       4                   4             198             1                       1                   2             200             0                       7                   7             201             0                       1                   1             202             1                       0                   1             203             1                       0                   1             206             0                       1                   1             209             0                       1                   1             213             0                       4                   4             220             0                       2                   2             221             0                       3                   3             225             0                       2                   2             227             1                       1                   2             229             0                       5                   5             231             3                       4                   7             232             1                       2                   3             234             0                       2                   2             235             0                       5                   5             237             0                       3                   3             242             0                       3                   3             245             0                       1                   1             255             0                       2                   2             260             0                       3                   3             264             0                       5                   5             266             0                       1                   1             269             1                       1                   2             274             1                       0                   1             282             1                       0                   1             286             0                       2                   2             287             0                       3                   3             295             0                       1                   1             298             0                       2                   2             300             1                       1                   2             336             0                       2                   2             360             0                       2                   2             364             0                       1                   1             368             1                       0                   1             370             1                       0                   1             375             1                       0                   1             393             0                       1                   1             415             0                       1                   1             428             0                       1                   1             432             0                       3                   3             452             0                       1                   1             454             0                       1                   1             485             1                       0                   1             496             0                       1                   1             500             1                       0                   1             506             0                       1                   1             534             1                       30                  31             567             1                       0                   1             572             0                       1                   1             574             0                       1                   1             605             0                       2                   2             607             0                       2                   2             684             0                       1                   1             698             0                       6                   6             708             0                       1                   1             714             0                       2                   2             727             0                       1                   1             732             1                       0                   1             764             0                       1                   1             809             0                       1                   1             828             0                       1                   1             837             0                       1                   1             883             0                       3                   3             970             0                       1                   1             999             1                       33                  34              **1) I'm quite positive I can't transform this into normalized data.**   I have performed, in order from most normalized to least normalized (according to my software): Box-Cox, Loglogistic (2P), Johnson (Type SU Unbounded), Lognormal (2P), Beta 4P, Weibull (2P), Gamma (2P), Weibull (3P), Logistic (2P), Largest Extreme Value (2P), Loglogistic (3P), Lognormal (2P), Normal, Gamma (3P), Exponential (1P), Exponential (2p), Half-Normal, (1P), and Half- Normal (1P) transformations. Regardless, the AD p-value in all of these cases was less than .05 and usually much, much less. From what I've learned, this is bad and means the transformation cannot be applied. I may as well give a shout out to the software add-in I used in Excel for being able to do all of these transformations for me at once and give me the results in graph and numeric form: SigmaXL - whom I'm not affiliated with at all nor do I know if they did all these correctly as I have not heard of 75% of them. I'm just assuming they know what they're doing as they seem credible and have pretty helpful & lengthy videos. The videos do leave alot of holes unfilled despite their usefulness though - as you may be able to tell from my question(s).   I also tried some manual transformations on my own like proportions involving subgroups or taking the ln of each value - these didn't help - or I did them incorrectly. **Any suggestions here? Am I completely wrong and I actually could transform this into normalized data?** **2) Non-parametric "transformation"?, "ranking"?, or what? And what do I do with those darn tests everyone speaks of. Testing is great, but what ACTION do I take?** Reading the above question #1, it seems apparent I'm going to have to do some sort of non-parametric "transformation" or "ranking." The problem is, as I stated earlier, all I can find are non-parametric "tests." I can't seem to find any way to transform the data using a non-parametric technique, nor do I know what to do with the values from some of the tests I've run like Mann Whitley. It will say the test is significant and then give me a Mann Whitley value, but I can't figure out what that means for instance, find out what to do with it through research on the net, and thus how to practically apply it. Am I even supposed to try to fit non-parametric transformed data (if its even possible to transform non-parametrically) with normalized parametric data into a binary logistic regression? **What do I do here?** **3) Blanks** As I stated before, I decided to, as part of my process, scrub out the blanks as they indicated 95% of the time the zip code was wrong anyway. I do realize this may screw up how I transform certain data columns, BUT I actually do want everything to fit into and be calibrated towards a blank free model anyway. I expect and hope for disagreements/suggestions against this. **Thoughts and Opinions Welcome** As a side note: No need to question the validity of the data I am comparing my returned data from x categorical input to. For the sake of complexity, just assume it is correct.


title: Is there any shortcut in SAS to interpret Limited Probability Model coefficients?

text: By LPM I mean that the dependent variable is polychotomous (e.g. 1,2,3 4) and (NOT binary 1 or 0). I know how to transform the coefficients manually by reverse calculating the PDF. Is there any command in SAS that would do it automatically? If I start let's say               proc logistic     OR         proc probit      I want to explain what happens if the independent variable changes by 1 unit?


title: Computing P-value

text: I am doing a lower-tail test and my z-statistic is -1.57. Now I have to calculate the p-value. My way of finding the p-value using a z-table is: I looked for row -1.5, then look for column 0.07. The entry is 0.0582. Since the test is a lower-tail form, the p-value is 0.0582. But from this source I found the p-value for the given problem is 0.58051. How was this calculated?


title: Difference between "probabilistic modeling" and "statistical modeling"

text: Is there a formal difference between the terms "Probabilistic Model" and "Statistical Model"? Is there a methodological difference between the two or just a preference of terminology? I see the former used most often with graphical models and Bayesian models, while "statistical" seems to be anything specifically Frequentist.


title: multivariate integration in closed form

text: I am looking for a closed form for the below integral but since I don't have the necessary backgrounds I am not able to solve it: i know the final solution is in the form of modified Bessel functions and polynomials and also I know it might be related to moments of multivariate inverse gaussian distributions: here is the integral assuming z and y are vectors. and $\lambda$ and $\Gamma$ are matrices. $\Gamma$ can be diagonal though. $U$ is a diagonal matrix containing $\sqrt z_i$ as its diagonal elements. $$ U = \begin{bmatrix} \sqrt{z_1} & \cdots & 0 \\\ \vdots & \ddots & \vdots \\\ 0 & \cdots & \sqrt{z_d} \end{bmatrix}$$ $$ \lambda = \begin{bmatrix} \lambda_{11} & \cdots & \lambda_{1d} \\\ \vdots & \ddots & \vdots \\\ \lambda_{d1} & \cdots & \lambda_{dd} \end{bmatrix}$$ $$ \Gamma= \begin{bmatrix} \Gamma_{1} & \cdots & 0 \\\ \vdots & \ddots & \vdots \\\ 0 & \cdots & \Gamma_{d} \end{bmatrix}$$ $$ \int^{\infty}_{0} \frac{1}{\sqrt{z_1 ... z_d}} \exp^{(-\frac{1}{2} y^T(U \Gamma U^T)^{-1} y)} \exp^{(-\sqrt{z}^T\lambda\sqrt{z})} \ \ \mathrm{d}z_1...\mathrm{d}z_d$$ thank you


title: Regression by subset in R

text: Perhaps this is more of a programming question than a stats question, but I'm sure someone here has the answer. I am new to R and am trying to run a linear regression on multiple subsets ("Cases") of data in a single file. I have 50 different cases, so I don't want to have to run 50 different regressions...be nice to automate this. I have found and experimented with the `ddply` method, but this, for some reason, returns the same coefficients to me for each case. Code I'm using is as follows: `ddply(MyData, "Case", function(x) coefficients(lm(Y~X1+X2+X3, MyData)))` Results I get, again, are the same coefficients for each "Case". Any ideas on how I can improve my code so that the regression runs once for each case and gives me unique coefficients for each case?


title: Performing comparisons after a cluster analysis

text: I have a data set of about 40 subjects, and 10 variables, some measuring performance on a task, and some measuring some characteristic. Each of the performance variables consists of the proportion of correct responses (range 0 to 1.0) to some set of questions. The characteristic variables (such as verbal IQ) are ordinal, and may or may not follow a normal distribution. The subjects have been bifurcated into two groups by way of cluster analysis. I'd like to compare the two clusters on each performance and characteristic measure to test if they are significantly different from each other on that measure. However, the clustering process entails a violation of assumptions that need to be met for standard statistical tests. What sort of tests should I use to test if the two clusters differ on the measures? Thank you.


title: Posterior in latent Dirichlet analysis

text: I have a question regarding LDA (Latent Dirichlet Analysis) - what is the correct formulation of the posterior? In http://www.cs.princeton.edu/~blei/papers/Blei2011.pdf it is $p(\beta_{1:K}, \theta_{1:D}, z_{1:D} | w_{1:D}) = \frac{p(\beta_{1:K}, \theta_{1:D}, z_{1:D}, w_{1:D})}{p(w_{1:D})}$, but in http://www.cs.princeton.edu/~blei/papers/BleiNgJordan2003.pdf and http://www.cs.princeton.edu/~blei/papers/BleiLafferty2009.pdf it looks different. I would be very grateful if someone could explain me where do these differences come from.


title: Analysis of repeated measures when fixed effect changes within-subjects?

text: I'm helping a colleague with a repeated-measures analysis in which each subject's fixed effect changes over the course of data collection. He's measuring a continuous, positive response variable (Response) in agricultural plots (ID) (which themselves are grouped within Blocks) as a function of Crop type. The response variable is measured 4 times over the course of the season. The analysis would normally look something like this:               modelout <- lmer(Response ~ Crop + (Month|ID) + (1|Block) , data=cropdata,                       family=gamma(link="inverse"))      However, just to make things complicated, the first measurement is done prior to planting, so it's influenced by the _previous year's_ crop type. That is, for plot $A_1$, Month 1's Crop is Wheat, but Month 2-4's Crop types are Barley. Is there any way to analyze this in a single model? Alternatively, I see two options:   1. Analyze Month=1 separately from the repeated-measures analysis for Month=2-4;    2. Use compound fixed effects, e.g., Wheat-Barley, Wheat-Corn, etc. Any suggestions? Edited to add: another thought - can we just add the first month's response and crop types to the model as covariates (continuous and categorical, respectively)?


title: How to understand this objective function in deep learning

text: I'm going through Christopher Manning's tutorial from NAACL 2013 "Deep Learning for NLP (without Magic)" and he gets to the point where he's showing how to do unsupervised pre-training. He's saying that we compute two different scores $s$ and $s_c$ and then we want to train the network so that $s > s_c$. He claims that the objective function to use is $$ J = \text{max}(0, 1 - s + s_c) $$ I don't see how minimizing this function ensures that $s > s_c$. It seems to me that the value is greater than zero whenever the sum of $s$ and $s_c$ is greater than 1 and it doesn't matter which one is greater than the other. Could someone help clear this up for me? What am I missing? The video is here: http://techtalks.tv/talks/deep-learning-for-nlp-without- magic-part-1/58414/ and the part I'm asking about is at about the 50:30 mark. The slides are at http://nlp.stanford.edu/courses/NAACL2013/NAACL2013-Socher- Manning-DeepLearning.pdf and the slide in question is number 50.


title: How do I train HMM's for classification?

text: So I understand that when you train HMM's for classification the standard approach is:   1. Separate your data sets into the data sets for each class   2. Train one HMM per class   3. On the test set compare the likelihood of each model to classify each window But how do I train the HMM on each class? Do I just concatenate the data pertaining to one class together? But isn't the time series data meant to be sequential - and if I do that then I am saying that some data points are consecutive when they are not? To be more concrete I have some EEG data which is a 96xT matrix where I have 96 feature vectors which are the power spectral densities of different frequencies from different channels and T is the length of time of the signal (at some sampling rate) This can be divided into windows which I know from the experimental protocol (the data is labelled) and so I can gather together sets of 96*t matrices for each class. Where t is less than T and denotes the size of each window. How do I then train the HMM on this data? If it helps I am trying to use the pmtk3 toolkit, but I open to using anything really - it just has to be able to deal with real-valued observations as the power spectral densities are continuous not discrete (the default MATLAB toolbox can only deal with discrete observations). The aim is to be able to classify EEG data windows to a given mental state having trained upon the labelled data. It is a brain-computer interface problem using the Berlin BCI Competition data.


title: Dimensionality reduction by predominant component analysis

text: Does any one know the concept of dimensionality reduction by predominant component analysis. ? This is introduced in one of the scientific papers. Please could any one explain how it works. Below is the link for the paper. The concept in described in section 4.3 http://www.eyemaginary.com/Rendering/TurnColorsGray.pdf


title: Distance between two random vectors

text: I have two random vectors, $A$ and $B$ with each consisting of $n$ geographical co-ordinates $(x_1,y_1),(x_2,y_2)\dots (x_n,y_n)$ and $(\tilde{x}_1,\tilde{y}_1),(\tilde{x}_2,\tilde{y}_2)\dots (\tilde{x}_n,\tilde{y}_n)$, respectively. Suppose we have a given distribution; let $p_1,p_2,\dots,p_n$ be the probabilities associated with points in $A$ and $\tilde{p}_1,\tilde{p}_2,\dots,\tilde{p}_n$ be the respective probabilities for points in $B$. My question is: What is the best way to find the distance between $A$ and $B$? I am looking for a distance measure which would give me zero value if co- ordinates in $A$ and $B$ are the same (appearing in any order) and it returns a large value when the co-ordinates in $A$ and $B$ are far apart. I considered using a covariance matrix but I don't get the intuition behind using that. In wikipedia, http://en.wikipedia.org/wiki/Statistical_distance, there are a lot of statistical measures listed. Can anybody help me with insights on how to choose a proper distance measure to suit the specifics of my problem?


title: Selecting from a bigger control group for 1:1 ratio

text: What I have is a treatment group consisting of 90 individuals. I have a pool of individuals eligible for my control group, which I want to be matched by age and gender in a 1:1 ratio. Is there any practical package or function in R to find such optimal control group?


title: Prior on a non identifiable parameter-MCMC integration

text: To introduce the problem I will explain the Projected normal distribution. Let $\mathbf{z}_i=(z_{i1},z_{i2})$ be a bivariate vector distributed as a bivariate normal with vector mean $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$. If we make the following transformation $x = \text{atan2}(z_{i2},z_{i1})$ (the the definition of the $\text{atan2}$ can be seen at http://en.wikipedia.org/wiki/Atan2), $x_i$ is said to be distributed as a projected normal variable with parameter $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$. Since $\text{atan2}(z_{i2},z_{i1})=\text{atan2}(cz_{i2},cz_{i1})$ for every $c>0$, the parameter of the projected normal distribution are not identifiable unless we use an identification constraint. Generally the variance of the second component of $\mathbf{z}$ is assumed to be $1$, and we let $\mathbf{V}$ be $\boldsymbol{\Sigma}$ with the constraint and with $\boldsymbol{\mu}^*$ the associated vector of $\boldsymbol{\mu}$, from $(\boldsymbol{\mu}, \boldsymbol{\Sigma})$ we can always compute $(\boldsymbol{\mu}^*, \mathbf{V})$. Consider now to have $n$ observations $x_i, i=1,..., n$ and we want to estimate the posterior distribution $\boldsymbol{\mu}^*, \mathbf{V}| x_1,...,x_n$. I have to choose prior distributions, for $\boldsymbol{\mu}^*$ I use a normal distribution, but I have some problems with the prior on $\mathbf{V}$. In this simple case I can parameterize $\mathbf{V}$ in terms of variance and correlation and put priors on the variance of the fist component and the correlation coefficient but in some cases $x_i$ can be multidimensional and then $\mathbf{V}$ can be, for example, a $10 \times 10$ matrix and a parametrization in term of variance and correlation cannot be feasible. This is my question/idea: I can estimate the following posterior $\boldsymbol{\mu}, \boldsymbol{\Sigma}| x_1,...,x_n$ and then map the posterior samples to the posterior samples of $\boldsymbol{\mu}^*, \mathbf{V}| x_1,...,x_n$. I am not sure that I can do this but I think that it is just an MCMC integration, something like: $$ f(\boldsymbol{\mu}^*, \mathbf{V}| x_1,...,x_n) = \int_{(\boldsymbol{\Sigma}, \boldsymbol{\mu})\rightarrow(\mathbf{V}, \boldsymbol{\mu}^*)} f(\boldsymbol{\mu}, \boldsymbol{\Sigma}| x_1,...,x_n) d(\boldsymbol{\Sigma}, \boldsymbol{\mu}) $$ where $f()$ is the pdf and the integral $\int_{(\boldsymbol{\Sigma}, \boldsymbol{\mu})\rightarrow(\mathbf{V}, \boldsymbol{\mu}^*)} d (\boldsymbol{\Sigma}, \boldsymbol{\mu})$ must be intended as the multiple integral over all the values of $\boldsymbol{\Sigma}$ and $\boldsymbol{\mu}$ that can be mapped into $\mathbf{V}$ and $\boldsymbol{\mu}^*$ (I am not sure that the integral is well written). Since $\boldsymbol{\Sigma}$ is a covariance matrix, am I right in thinking I can use a Normal Inverse Wishart distribution on $\boldsymbol{\mu}, \boldsymbol{\Sigma}$ and the problem is solved? Moreover, how can i interpret the parameters of the Normal Inverse Wishart distribution


title: Plotting confidence intervals for the predicted probabilities from a logistic regression

text: R and Statistics newbie here. Ok, I have a logistic regression and have used the predict function to develop a probability curve based on my estimates.               ## LOGIT MODEL:     library(car)     mod1 = glm(factor(won) ~ as.numeric(bid), data=mydat, family=binomial(link="logit"))          ## PROBABILITY CURVE:     all.x <- expand.grid(won=unique(won), bid=unique(bid))     y.hat.new <- predict(mod1, newdata=all.x, type="response")     plot(bid<-000:1000,predict(mod1,newdata=data.frame(bid<-c(000:1000)),type="response"), lwd=5, col="blue", type="l")      This is great but I'm curious about plotting the confidence intervals for the probabilities. I've tried plot.ci() but had no luck. Can anyone point me to some ways to get this done, preferably with the car package or base R. Thanks.


title: Inequality distribution function

text: > **Possible Duplicate:**   >  Inequality distribution function I am trying to understand a proof (no homework). If $\frac{\overline{F}(x-y)}{\overline{F}(x)} \rightarrow 1$ for $x \rightarrow \infty$, then why does the following follow: $\int\limits_{0}^{s}\frac{\overline{F}(t-u)}{\overline{F}(t)}dF(u)\geq 1-\frac{\epsilon}{4}$, where $s,t>0$.


title: About the GMRFlib

text: I am trying to access the GMRFlib as mentioned in the link http://www.math.ntnu.no/~hrue/GMRFsim/doc/html/. However, I think it is not available anymore. I couldn't see the tarball mentioned. Does anyone have access to it?


title: Simulating from a normal with "unknown" variance

text: Suppose I want to performing sampling from a normal distribution with an unspecified variance, and I want a way to sample so that I am in some sense "averaging out the possible values of the variance." In other words, I want to use the simulation results elsewhere, and hear results are highly dependent on the variance, and so I want to avoid "favoring" one choice of variance over another. I was thinking maybe I could put a prior on the variance and derive the marginal posterior of the variance given the draws, find the MAP estimate of the posterior, say, $m$, and then simulate again from a $\mathcal{N}(0,m)$ distribution? If not, how else could I do this?


title: sampling from confidence intervals

text: If we are provided a mean and (95%) confidence interval, is it possible to set up a system in which we draw random values that are outside the CIs 5% of the time? My _intuition_ is that if one can generate a distribution in which 68% of the random draws fall within a given interval, then it should also be possible to generate a distribution in which 95% of the draws fall within an interval. Examples in R always helpful!


title: Cross-validation in multi-level model

text: Suppose I want to estimate the out-of-sample prediction error of a boosted regression model that has random intercepts and slops. There are $G$ groups and $N$ observations. If I want to estimate the out-of-sample prediction error using $k$-fold cross-validation, how do I set up the data partitioning? Is it more complicated that $k$-fold cross-validation? Note: my use case here is the prediction of the data from a _new group_.


title: Hard Case - prediction of chain stores revenue

text: * Data about average monthly revenue from 2000 stores around whole country.   * Gini coeff. of reve around 20%, with 50% of observation around average, very thin tails of distribution   * Explanatory variables: geolocation of competitors, traffic drivers, purchasing power, population etc. in cachment area (most scatter plots cloud-shaped with small correlation)   * some outlier cases I've tried glm, svr, regression trees with boosting. All models predictions are only a little bit better then average (RMSE in avg model = 40k, RMSE in regression = 35k, reve avg = 130k), with small errors in centre and huge overestimation in left tail and underestimation in right tail. Any idea for new model type or prediction technique that can fit better in tails without big disturbance in the center ?


title: conditional independence in repeated measures design

text: How the responses are independent when conditioned on random effect in repeated measure analysis (linear mixed model)?


title: Software or libraries to create doc-term matrix

text: does anyone know some **_Java libraries_** to create the document-term matrix for a large number (50,000) of documents ? I wish this library encompasses preprocessing functionalities, like stop-word and punctuation removal, stemming, etc. What's more, I wish to use TF*IDF weighting scheme, and normalization functionality. Preferably, I prefer Java libraries for the convenience of development. Thanks very much for any recommendation.


title: non parametric index model (Klein Spady) for more than two choices in R

text: The np package offers an easy way to fit the Klein Spady estimator for index model on dichotomous responses. How to use the np package to extend these methods to the polychotomous (i.e., more than two possible responses) case, e.g., Lee, journal of econometrics 1995?


title: Making boxplots of hourly data in R

text: I have air temperature measurements from two fixed locations measured at hourly intervals. The code below is a random set of numbers used to represent the format of my data:               set.seed(1)     RandData <- rnorm(8760*2,sd=10)     Locations <- rep(c('UK','France'),each=8760)          Date = seq(from=as.POSIXct("1991-01-01 00:00"),                    to=as.POSIXct("1991-12-31 23:00"), length=8760)          Final <- data.frame(Loc = Locations,                         Doy = as.numeric(format(Date,format = "%j")),                         Tod = as.numeric(format(Date,format = "%H")),                         Temp = RandData)      I can plot the variation in temperature as a funtion of day of year with the following code:               require(lattice)     xyplot(Temp~Doy | Loc, data = Final, col = "black", type = "l")      This would show the annual pattern of the data. However, what I would like to do is to produce boxplots of the variation in temperature for different times of the day. So, for the example above I would like two figures, one for each country and each figure should be composed of box plots showing the variation in temperature at 00:00, 01:00... and so on, referring to Final$Tod. How can this be achieved? Many thanks for your help.


title: Is there a quantitative way to compare the distribution shape of different samples?

text: I am conducting some research which involves visually/graphically observing the differences between the shapes of the distributions of different samples. I would like to automate this process (at least somewhat), so that I can scale the number of samples I look at (as well as speeding things up, reducing human error etc.). Is there a way to _quantitatively describe/measure the shape of a distribution_ so that comparisons between shapes can be made algorithmically?


title: Interpretation of multiple logistic regression with interactions in R

text: I am trying to look at whether 2 variables (one dichotomous categorical and one continuous) predict the occurrence of a dichotomous categorical dependent variable.               dependent variable is LENIpos - 0 = no event, 1 = event     predictor variables are Hip.Prox.Femur - 0 = no hip fracture, 1 = hip fracture                     and     age (continuous)      Both predictor variables have significant p values in separate chi square test and Mann Whitney U test respectively. When I run a logistic regression `glm(LENIpos ~ age + Hip.Prox.Femur, family = "binomial)`, the variables come out as not significant. (1) However, when I run the logistic regression with interactions `glm(LENIpos ~ age * Hip.Prox.Femur...)` (2), they are no both significant. How is this to be interpreted? Example R outputs: (1)               Call: glm(formula = LENIpos ~ age + Hip.Prox.Fem, family = "binomial",          data = dvt)          Deviance Residuals:           Min       1Q   Median       3Q      Max        -0.9346  -0.7826  -0.4952  -0.3374   2.1897            Coefficients:                              Estimate Std. Error z value Pr(>|z|)        (Intercept)              -3.46888    1.00693  -3.445 0.000571 ***     age                       0.02122    0.01519   1.397 0.162535       Hip.Prox.Femhip fracture  0.72410    0.57790   1.253 0.210212              (Dispersion parameter for binomial family taken to be 1)              Null deviance: 145.23  on 151  degrees of freedom     Residual deviance: 135.48  on 149  degrees of freedom     AIC: 141.48          Number of Fisher Scoring iterations: 5      (2)               glm(formula = LENIpos ~ age * Hip.Prox.Fem, family = "binomial",          data = dvt)          Deviance Residuals:              Min       1Q   Median       3Q      Max           -1.0364  -0.7815  -0.5373  -0.1761   2.3443            Coefficients:                                  Estimate Std. Error z value Pr(>|z|)       (Intercept)                  -5.89984    1.98289  -2.975  0.00293 **     age                           0.05851    0.02818   2.076  0.03788 *      Hip.Prox.Femhip fracture      5.04990    2.46269   2.051  0.04031 *      age:Hip.Prox.Femhip fracture -0.06058    0.03339  -1.814  0.06965 .                (Dispersion parameter for binomial family taken to be 1)              Null deviance: 145.23  on 151  degrees of freedom     Residual deviance: 131.82  on 148  degrees of freedom     AIC: 139.82          Number of Fisher Scoring iterations: 6


title: What is this problem called in machine learning research?

text: > Given 2 sets of entities from 2 different classes described by properties > $f_1...f_n$, any 2 entities from the 2 different classes together have a > score: ${\rm Score}(c_1e_1,c_2e_1) = y$. How can one derive a function / > model predicting a score of any pair of the two different classes given its > properties? What's the term for this problem in research?


title: Data on a circle - testing if points are close to another given point

text: I have a set of points that lie on a circle, and I know their angles. Is there a way that I can test whether they are closer to another point `x` on the circle than expected by chance (assuming a uniform distribution on the circle)? Can I calculate mean arc length between `x`and the data, and then calculate the same for each of many simulated data sets generated by the uniform distribution?


title: Limited dependent variable → ordered probit model with Stata

text: I have a dependent variable, credit rating, and it is a limited dependent variable. Independent variables are profitability, firm size, stock return, etc. I was advised to use Stata and run regression with fixed effects and then `predict yhat` in order to obtain a target rating that a certain firm wants to target in this way: `xtreg dependent independent, fe` with regression equation below: Credit rating = $\alpha + \beta_1$R&D $+\beta_2$Size $+\beta_3$Selling Expense $+\beta_4$Market to Book $+\beta_5$Profitability $+\beta_6$Operating Risk $+\beta_7$Tangibility + error term Then run `predict yhat` in order to get a proxy for the credit rating a certain firm wants to target. The dependent variable credit rating is limited, so is the panel regression for the model: `xtreg, fe` oK, or should an ordered probit model be used instead? (In other studies with credit ratings they use ordered probit.) And if a probit model should be used, what command is the corresponding command `predict yhat` (in the case of `xtreg fe`) for the ordered probit model? Or in other words, how do I get the rating that a firm wants to target (from the panel regression with `xtreg, fe` and then `predict yhat` for the fitted values) from the ordered probit model?


title: How to deal with bouts of equilibrium in an online learning setting?

text: I have a kalman filter (a recursive least square filter, really) regressing over real-time streams of data. Because the data-generating process varies slightly over time I add an exponential forgetting factor to the kalman filter.   It all fits nicely with good prediction. The problem I have is that the variables in input streams sometimes enter long periods of equilibrium. In those periods of time all I receive in the input streams are the same $y_t$ and same $x_t$ (except for background white noise). If the equilibrium remains for long enough the kalman fit starts degrading as it forgets old useful data and gives more and more importance to the background noise in equilibrium. What I'd like to know is if there is any method to avoid overfitting the filter to useless new data. Right now I am using very ad-hoc methods (tallying a moving standard deviation to use as weight in the regression), but there must be better way to regularize a regression in an online setting.


title: Learning models for dependent point processes

text: This is largely a literature request. I have some data consisting of, for each use of a credit card, the card owner, the store, and the time of the use. This could be viewed as a collection of dependent point processes associated with the edges of a big graph, where the vertices are credit card owners and stores that accept credit cards. I'm interested in trying to model this information, in particular in such a way that one could try to simulate data from the model afterwards without too much difficulty. It is ok if one can only simulate at small discrete time steps. Is there a literature on this type of topic? Any great papers to read?


title: What is the influence of C in SVMs with linear kernel?

text: I'm currently using an SVM with a linear kernel to classify my data. There is no error on the training set. I tried several values for the parameter C (10^-5, ..., 10^2). This did not change the error on the test set. Now I wonder: is this an error caused by the ruby bindings for libsvm I am using (https://github.com/febeling/rb-libsvm) or is this theoretically explainable? Should the parameter C always change the performance of the classifier?


title: Bag of words vs vector space model?

text: What is/are the difference/s between these text representation models: Bag of words and vector space model?


title: Simulating from Posterior Predictive Over Many Periods

text: Suppose I obtain a posterior distribution through MCMC. This will give some simulated values for the parameters to a statistical model. If I calculate the posterior predictive distribution one period ahead, then I can simply use the simulated posterior distribution values to make the prediction. However, if I want the posterior predictive distribution in each of $t$ periods ahead, then I have some options. I could:   1. Use the parameters simulated from the posterior distribution to simulate from the posterior predictive in each of the $t$ periods. That is, you could simulate from each simulated posterior over the $t$ periods (saving each period).   2. Simulate a different posterior for each $t$ instead of only once for each $t$. Use the period-specific posterior parameters to generate the posterior predictive.   3. Scale the posterior to what would be appropriate for each period and then simulate each period directly. For instance, if the relevant parameters are the mean and variance of a normal distribution, then they can be multiplied by $t$ to obtain the posterior distribution at each $t$. This direct approach works well if I only care about the posterior predictive at time $t$, but may not generally be as effective if I care about the intervening periods. I have been taking the first approach, but I'm not sure if/why it is better than the second approach.


title: mean difference for count data

text: I have two samples s1 and s2 of count data. The sample size is > 1000 each. The distributions look similar to a Poisson distribution but the variance is much larger than the mean. How do I test whether the mean of s1 is larger than the mean of s2?


title: What are some valuable Statistical Analysis open source projects?

text: What are some valuable Statistical Analysis open source projects available right now? Edit: as pointed out by Sharpie, valuable could mean helping you get things done faster or more cheaply.


title: Unbiased hypothesis tests

text: * Is there some textbook or expository account showing that the definition of "unbiased test" bears the same sort of relation to "unbiased estimator" that interval estimation generally bears to confidence intervals? (One rejects $H_0 : \mu=\text{whatever}$ precisely if $\text{whatever}$ is not within a confidence interval for $\mu$.)   * Similarly, an account of the relationship between the concept of "uniformly most powerful test" and that of "uniformly most accurate confidence interval"?   * If I'm not mistaken, the F-test of $\sigma_1^2=\sigma_2^2$ versus "$\ne$", where these are variances of normally distributed populations, is not unbiased. Does that actually make it an inferior test, and what improvements on it are thereby called for?


title: How high must logistic covariates' predictive accuracy be for a reversal effect to show up?

text: I am modeling an outcome for hospital patients, 'RA' (whether readmitted; 0=No, 1=Yes). My predictor of interest is 'HHS' (whether referred to Home Health Services such as from a visiting nurse; 0=No, 1=Yes). Those referred readmit at a 15.2% rate; others, 9.2%, but the former are needier, sicker patients. Conventional thinking is that if we controlled for severity of illness this difference would not only be washed out but would reverse itself. In other words, holding constant the severity of illness, having HHS should mean a lower RA rate. With HHS as the sole predictor, its coefficient (B) in a **logistic** regression = 0.6 (N ~ 25k). B is reduced to 0.2 with a group of covariates controlled, each accounting for some aspect of severity of illness, but B doesn't fall below zero. HHS alone explains only about 1% of the variance in RA; with the other predictors, this becomes 4%.* Perhaps this is the problem--that these covariates are not explaining enough variance to "succeed" in reversing the sign of the coefficient of interest. If this is true, is there a way to estimate how high **their** explained variance needs to be for such a reversal to show up? * * * *Using either of 2 pseudo-R-squared formulas; Cox & Snell's or Menard's [-2LL0 - (-2LL1)] / [-2LL0.]


title: Why are directed graphical models called Bayesian?

text: Why are directed graphical models called Bayesian,while those undirected not? Does that mean that Bayesian analysis is involved in the directed ones, and not in those undirected?


title: Determining variance of an U.M.V.U.E

text: Given $X_1...X_n$ i.i.d. Bernoulli R.V.s with parameter $\theta$, I found an UMVUE for $\tau(\theta)=\theta^2+\theta$ with Rao-Blackwellization. That process seemed fairly straight-forward to me, and I ended up with $\frac{T(n+T-2)}{n(n-1)}$, where T is the number of successes and n is the number of trials. Second, I calculated the Cramer-Rao lower bound and wound up with $\frac{(4\theta^2+4\theta+1)(\theta(1-\theta))}{n}$. I'd like to plot these both and see if the UMVUE attains the CRLB, but to do that I need to find the variance of my UMVUE! I imagine I want to go about this as $E([\mathrm{umvue}]^2) - E([\mathrm{umvue}])^2$, and I already know that $E([\mathrm{umvue}])$ is $\theta^2 + \theta$, since its unbiased. However the second moment is confusing me greatly; in fact I have no idea how to approach it. Am I missing something obvious? Thanks for the help! (And please let me know if there is anything I should clarify!) I finished up the question, and it did indeed involve the 4th and 3rd moments of the binomial. Should anyone ever encounter a similar problem, http://mathworld.wolfram.com/BinomialDistribution.html contains the 1st through 4th moments of the distribution, which should save you a great deal of computation!


title: Ordinal survey data or not?

text: We wish do some statistical analysis on the relation between being an entrepreneur and drinking and are wondering if we can consider the data derived from a survey as ordinal in nature. And if we can rank them as such. The possible answers are as follows:   * No, I have never drunk alcohol   * No, I no longer drink   * Yes, but only rarely   * Yes, I drink alcohol 2 or 3 days per month   * Yes, I drink alcohol 1 or 2 days per week   * Yes, I drink alcohol 3 or 4 days per week   * Yes, I drink alcohol 5 or 6 days per week   * Yes, I drink alcohol everyday The yes answers are in our opinion ordinal in nature. The problem is with the 2 no answers. Is it possible to rank the answers as we did here, or is it better to just group them as non-drinkers? EDIT: Research question: Do entrepreneurs drink more frequently than non- entrepreneurs (wageworkers)?


title: Determining whether a website is active using daily visits

text: ### Context: I have a group of websites where I record the number of visits on a daily basis:               W0 = { 30, 34, 28, 30, 16, 13, 8, 4, 0, 5, 2, 2, 1, 2, .. }      W1 = { 1, 3, 21, 12, 10, 20, 15, 43, 22, 25, .. }     W2 = { 0, 0, 4, 2, 2, 5, 3, 30, 50, 30, 30, 25, 40, .. }      ...     Wn       ### General Question:   * How do I determine which sites are the most active? By this I mean receiving more visits or having a sudden increase in visits during the last few days. For illustration purposes, in the small example above W0 would be initially popular but is starting to show abandoning, W1 is showing a steady popularity (with some isolated peak), and W3 an important raise after a quiet start). ### Initial thoughts: I found this thread on SO where a simple formula is described: ` // pageviews for most recent day   y2 = pageviews[-1]   // pageviews for previous day   y1 = pageviews[-2]   // Simple baseline trend algorithm   slope = y2 - y1   trend = slope * log(1.0 +int(total_pageviews))   error = 1.0/sqrt(int(total_pageviews))   return trend, error   ` This looks good and easy enough, but I'm having a problem with it. The calculation is based on slopes. This is fine and is one of the features I'm interested in, but IMHO it has problems for non-monotonic series. Imagine that during some days we have a constant number of visits (so the slope = 0), then the above trend would be zero. ### Questions:   * How do I handle both cases (monotonic increase/decrease) and large number of hits?    * Should I use separate formulas?


title: How should I distribute a classifier to customers?

text: When consulting, I often do my exploratory analysis and prototyping in R, and deliver results on the initial dataset to the client. The client wants to use the trained classifier in a production setting. For example, assume that the client wants to run 1000 records per day through the classifier. **How should I deliver the trained classifier?** For example:   1. Send paper describing the results and R code that will exactly replicate my work and hope that the client will be able to make a production server and implement it properly?   2. Abandon R after the proof of concept, and write up the classifier in a different programming language (e.g. python) that is more readily deployable by the client?   3. Try to set up a webserver that will run an instance of R?   4. Accept data from the client and run it upon request?


title: Cross validation using experimental design methods?

text: I have a question about using experimental design methods for cross validation in a non-conventional approach. I will use an illustration to explain the question. Imagine you have 100 time series stacked on top of each other and segmented into adjacent interval partitions (we can use 1990, 1991, 1992 for example). Then suppose we take the time series from the first year and transform them into some factor with levels. So, let's say A, B, C represent top, middle, and bottom 33% of performers, respectively. Then, we look at a continuous response variable such as percent returns over the next time interval. Ignoring any limitations for the moment, like unequal variance, we could use ANOVA to tell us whether or not a set of response variables have any statistically different means. Further, imagine we find the response variable corresponding to B is different than A and C. Now imagine we slide the box forward, and we apply the same method to the next time interval. I would like to have a test that looks at whether the properties of the response variables have the same statistical differences as the prior interval ... And possibly whether the underlying time series (block?) properties have changed. Imagine a two factor interaction plot with two overlapping V formations representing zero statistical difference. What is the appropriate approach to answer this, is it a variation of ANOVA (e.g. 2-Way)? More importantly, is there a specialized area dedicated to looking at the mixed signal time series cross validation perspective from more of an experimental design perspective? Somewhere along the lines of longitudinal analysis combined with cross-sectional analysis, perhaps?


title: Could somebody explain to me what this ARIMA model output says?

text: I've asked a few questions here before regarding my thesis. Although I try my best to follow-up on your suggestions, my statistical knowledge is limited but I try my utmost. Adding a predictive model to your thesis is not required (since it isn't taught during the studies) but my thesis coach insists. So I've just let SPSS dictate the best-fitting ARIMA model for my thesis. Basically, I have taken some internet data (`hbVol0LN` is number of tweets, `hbBullQuality0` is the ratio for postive against negative tweets, etc.) for 100 companies over 103 days. Here, the dependent variable is the return of the stock of each of those 100 companies per day. I already performed an OLS (although it has been pointed out that this is not the ideal model for my research, it is accepted by my coach), but now I believe this ARIMA model should hold the predictive value of the data. It is very hard to find annotated ARIMA output online, or a paper which describes the output in a way I can understand. Could you perhaps give me some insights of what this output is telling me? Any help at all is greatly appreciated. If you're having a hard time reading the graph, here's the full-size one: http://i.stack.imgur.com/H42YP.png Again, I cannot express how frustrating it is for somebody who has hardly had to do any statistics during his studies, having to produce a predictive model. Therefore, really, any help is appreciated. ![enter image description here](http://i.stack.imgur.com/H42YP.png)


title: Is it ok to clean up data from different conditions independently?

text: I have a dataset, which contains measurments from many different conditions. Since my hypothesis suggested a large difference for the measurements in each condition, in order to clean the data, I analyzed all conditions independently. I.e. I grouped all measurments from one conditions into quartiles, calculated the interquartile range and then removed any datapoint, which was further than 1.5 times the size of the interquartile range from the median. Now someone looked over the numbers and remarked, that the rate of removed data was much higher than usual. I re-checked my calculations several times, and arrived at such a high ratio of outliers each time. Now I was thinking, that if I had analyzed all measurements from each condition together, most likely much less measurements would have been removed, due to the differences between the conditions. However as far as I understand these cleaning methods, they are meant to be used on single (gaußian) distributions, and not the sum of two or more distributions. So which method for data cleaning would actually be the correct one: Cleaning all measurements together or cleaning each condition separately?


title: What is the expectation of exponential of the product of two random variables?

text: I am looking for examples of probability distributions that would allow me to characterize the distribution (at least approximately) and to compute the first two moments exactly of: $$ e^{aXY} $$ where $a$ is a positive real and $X$ and $Y$ are scalar random variables defined on the positive real line and jointly distributed with mean/variance $(\mu_{X},\sigma^2_{X})$ and $(\mu_{Y},\sigma^2_{Y})$ and covariance $\sigma_{XY}$. Any suggestions very welcome. I'm looking for examples, so hopefully there should be more than one approach to this.


title: Minimal state vector for a Kalman filter

text: I'm trying to work my way through a text on robotics and specifically trying to make sense of Kalman filters. I've used them before but I'd like to be able to understand and make my own for new systems. The question I have asks for the minimal state vector for a car with position X and velocity V at time t. The acceleration (A) is assigned randomly at each point in time. I have to come up with a minimal state vector so the the resulting system is Markovian. Here's my answer so far: The initial state vector should just be S = [X,V,A] since those are the three components that were given? The Markovian part is confusing to me. This means that the state vector can only have components that do not depend on past measurements? In which case it should only be S = [X,A] since acceleration was given to be random and position has no time dependence?


title: How can I cluster data in a grid-like fashion and heat map the averages in R?

text: I have a data frame of 3 columns. The first one is the response variable the second and the third ones are some criteria. You can create your own example similar to mine, using this piece of code with one difference; I have 120k of these rows.                   n<-10             data.frame(response=runif(n),x1=round(runif(n,min=0.2,max=3.8),2),x2=round(runif(n,min=14,max=180)))                response   x1  x2     1  0.007240072 0.99  94     2  0.585625664 3.26 175     3  0.060195378 1.52 153     4  0.806096047 1.90  15     5  0.715590971 2.87 161     6  0.840640566 3.06  73     7  0.757785139 3.38 125     8  0.835112330 1.43 158     9  0.588479082 1.68  59     10 0.963268147 0.54 108      **I would like to cluster the response column using x1 and x2 as the sides of a grid (non overlaping rectangular clusters while covering the whole population).** So any cluster will have 4 boundaries (min x1 - max x1 - min x2 - max x2 for each cluster) or just the vector of split points (x1.1, x1.2, x1.3 and x2.1, x2.2, x2.3 for nine rectangular cluster). Just like a basic Mondrian painting ![mondrian](http://i.stack.imgur.com/CbBLw.jpg) **And I would like to plot a heat map using the averages of the response column in each cluster or the number of instances in each cluster.** For the first part I used ctree function from party package. My attempts to use rpart function from the rpart package is thwarted by the 'intimidating' size of the data. See my question about this problem I managed to get the terminal nodes but not the node conditions. And how do I heat map such an object or data if I am given the split values?


title: Interactions with R? ANCOVA in R? I want to see if certain variables affect a dependent variable in R using repeated-measures

text: I have done an experiment in which the same people have been measured in two conditions, the same room with and without plants. I know they performed significantly better on a reverse digit span task in plant condition (it was relatively normally distributed data so I used a paired t-test), but want to see if the increase is related to mood (as measured by a likert-type scale) and their responses on a measure that assesses the "restorativeness" of the environment (another likert-type scale). Theoretically one should perform better on a digit span in a more "restorative" room. I might want to test if assessment of environment affects mood, which therefore affects performance on task. I would want to somehow compare this interaction (if that's the right term) in the plant condition and the no-plant condition. If that's too complex maybe I'd just want to see if the assessment of environment is a strong factor in their performance on the digit span. Is ANCOVA the right test for this? If so how should it be set up and if you'd like to be really generous how would I code this in R? Thanks for your help.


title: Estimating the spectral density

text: I'm interested in using the spectral density to determine the period of a time series $x_t$. According to Wikipedia, we can use the periodogram to estimate the spectral density. Going down the chain, Wikipedia says the periodogram is implemented using the FFT, but the details are not given. Is it literally the FFT on $x_t$? **Edit:** In R, there is a function spec.pgram that returns the raw periodogram. I tried to compare this to the FFT of the ACF, but I'm still not getting matching results (eventually I need to do this in Java, so I am trying to understand the implementation fully).               x = rep(c(5,1,5,10),4)     x.acf = acf(x)$acf[,,1]     Mod(fft(x.acf))^2           [1]  0.3906250  0.4757183  0.9898323 15.5956432  2.5543312  0.3962321      [7]  0.2068152  0.2068152  0.3962321  2.5543312 15.5956432  0.9898323     [13]  0.4757183      vs               x.pgram = spec.pgram(x, detrend=F)     x.pgram$freq          [1] 0.0625 0.1250 0.1875 0.2500 0.3125 0.3750 0.4375 0.5000          x.pgram$spec          [1]  3.8819277  3.4948335  2.9155061 69.0892857  1.5487796  0.9694522  0.5823580     [8]  0.1607143


title: Interpreting Granger causality test's results – i.e. $X = f(Y)$

text: The following is generated using http://www.wessa.net/rwasp_grangercausality.wasp:               Summary of computational transaction     Raw Input   view raw input (R code)     Raw Output  view raw output of R engine     Computing time  2 seconds     R Server    'Herman Ole Andreas Wold' @ wold.wessa.net          Granger Causality Test: Y = f(X)     Model   Res.DF  Diff. DF    F   p-value     Complete model  356              Reduced model   357 -1  17.9144959720894    2.94360540545316e-05          Granger Causality Test: X = f(Y)     Model   Res.DF  Diff. DF    F   p-value     Complete model  356              Reduced model   357 -1  0.0929541667364279  0.760632773377753      My interpretation is that: (i) Granger Causality Test: _Y = f(X)_ p-value = 2.94360540545316e-05   The p-value is very small, thus the null hypothesis _Y = f(X)_ , _X_ Granger causes _Y_ , is rejected. (ii) Granger Causality Test: _X = f(Y)_ p-value = 0.760632773377753   The p-value is near to 1 (i.e. 76%), therefore the null hypothesis _X = f(Y)_ , _Y_ Granger causes _X_ , cannot be rejected. Is my interpretation correct? [Related to full question in Interpreting Granger causality test's results.]


title: Checking model quality in linear regression

text: I found that in linear regression to check the model quality you can look at the plots described below (my questions are in bold).   1. scatter plot: plot Y against each X separately    2. scatter plot: plot residuals against each X in the model separately      * **Should it be scattered around the plot equally without any particular shape or form?**   3. scatter plot: plot residuals against each X not in the model separately (if linear include!)   4. scatter plot: plot residuals against Y and Yˆ      * **Should the datapoints be evenly scattered around the horizontal 0 line?**   5. normal QQ-plot of the residuals (check normality assumption.) Also: what does a plot of residuals vs the Y variable show and how should it look? Mine looks like this, does that seem ok? ![residuals vs y](http://s4.postimg.org/t3d5u394t/p1_residualsvsy.png) Grtz Stefan


title: Doubly robust estimation of causal effects implementation

text: Does anyone know of an implementation (other than the SAS macro) of the doubly robust estimation method found in: > Funk, M.J., Westreich, D. et al (2011). Doubly robust estimation of causal > effects. _American Journal of Epidemiology, 173_ (7): 761-767.[DOI]?


title: Time-series classification - very poor results

text: I am working on a time series classification problem where the input is time series voice usage data (in seconds) for the first 21 days of a cell phone account. The corresponding target variable is whether or not that account cancelled in the 35-45 day range. So it is a binary classification problem. I am getting very poor results from all of the methods that I have tried so far (to varying degrees). First I tried k-NN classification (with various modifications) and got extremely bad results. This lead me to extract features from the time series - i.e. mean, variance, max, min, total zeros days, total trailing zero days, difference between first half average and second half average, etc. and the most predictive features seemed to be total zeros days and total trailing zero days (using several classification algorithms). This performed the best but the performance was still not very good. My next strategy was to oversample the negative instances in my training set since there were so few of them. This resulted in more correct cancellation prediction but at the expense of more false-positives. I'm starting to think that perhaps the time series usage data itself is simply not very predictive (though common sense says that it should be). Perhaps there is some latent variable that I am not considering. Looking at the data also shows some strange behaviour. i.e. some examples show very little or decreasing usage (or sometimes none at all) and do not cancel, and some show a ramp up in usage that do cancel. Perhaps this contradictory behaviour does not generate a very clear decision boundary for a classifier. Another possible source for error is the fact that many training examples are very sparse (i.e. many days with 0 usage). One idea that I have not tried yet is to split the time series into segments and generate some features that way, but I do not have high hopes.


title: Fitting polytomous logistic regression with R

text: I'm running into troubles fitting a polytomous logistic regression model using grouped data. The data are of the form (dput at bottom):               > head(alligator)         lake  sex  size    food count     1 Hancock male small    fish     7     2 Hancock male small  invert     1     3 Hancock male small reptile     0     4 Hancock male small    bird     0     5 Hancock male small   other     5     6 Hancock male large    fish     4      And I've tried to fit the model with `vglm()` from package VGAM:               > result <- vglm(food~lake+size+sex, data=alligator, fam=multinomial, weights=count)     Error in if (max(abs(ycounts - round(ycounts))) > smallno) warning("converting 'ycounts' to integer in @loglikelihood") :        missing value where TRUE/FALSE needed     In addition: Warning messages:     1: In checkwz(wz, M = M, trace = trace, wzepsilon = control$wzepsilon) :       96 elements replaced by 1.819e-12      It was also suggested to look at `mlogit()` from package `globaltest` (on Bioconductor), but it does not appear to support grouped data. It obviously doesn't support the `weights` parameter, but I can't find where the equivalent parameter is documented:               source("http://bioconductor.org/biocLite.R")     biocLite("globaltest")          result <- mlogit(food~lake+size+sex, weights=count, data=alligator)     Error in mlogit(food ~ lake + size + sex, weights = count, data = alligator) :        unused argument(s) (weights = count)      If anyone could put me down the right path, I'd appreciate it!               > dput(alligator)     structure(list(lake = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L,      3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,      3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,      4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("George", "Hancock",      "Oklawaha", "Trafford"), class = "factor"), sex = structure(c(2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("female",      "male"), class = "factor"), size = structure(c(2L, 2L, 2L, 2L,      2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,      2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L,      1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,      2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,      1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L), .Label = c("large",      "small"), class = "factor"), food = structure(c(2L, 3L, 5L, 1L,      4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L,      2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L,      3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L,      5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L,      1L, 4L, 2L, 3L, 5L, 1L, 4L, 2L, 3L, 5L, 1L, 4L), .Label = c("bird",      "fish", "invert", "other", "reptile"), class = "factor"), count = c(7L,      1L, 0L, 0L, 5L, 4L, 0L, 0L, 1L, 2L, 16L, 3L, 2L, 2L, 3L, 3L,      0L, 1L, 2L, 3L, 2L, 2L, 0L, 0L, 1L, 13L, 7L, 6L, 0L, 0L, 3L,      9L, 1L, 0L, 2L, 0L, 1L, 0L, 1L, 0L, 3L, 7L, 1L, 0L, 1L, 8L, 6L,      6L, 3L, 5L, 2L, 4L, 1L, 1L, 4L, 0L, 1L, 0L, 0L, 0L, 13L, 10L,      0L, 2L, 2L, 9L, 0L, 0L, 1L, 2L, 3L, 9L, 1L, 0L, 1L, 8L, 1L, 0L,      0L, 1L)), .Names = c("lake", "sex", "size", "food", "count"), class = "data.frame", row.names = c(NA,      -80L))


title: Two R packages for topic modeling, LDA and topicmodels?

text: It seems that there have two R packages for running Latent Dirichlet Allocation. One is LDA, authored by Jonathan Chang; and another is authored by Bettina Grün and Kurt Hornik. What are the differences between these two packages, in terms of performance, implementation details and extensibility?


title: Compositional analysis vs. Manly selection ratios

text: I am comparing habitat availability vs. use by a bird species using the R `compana()` function (package adehabitat) which uses Aebischer et al. (1993) method. It is a 3rd order habitat selection analysis, so the availability of habitats is different for each bird. The `compana` randomisation test resulted in a nice habitat ranking plus a table with the significance of paiwise comparisons. The only thing I need more, are log-ratio differences or Manly selection ratio's (considering that these are given in almost every paper on habitat selection). Thus, I tried the `Wi` function that computes Manly's selection ratios. The problem is that the resulting ranking of habitats is very different from the one produced by `compana`, and the `compana` result seems to fit much better to the raw data. Apparently the results of these methods cannot be combined. So now I am a bit lost. I hope someone is willing to explain   * Why these two methods (compana and Manly selection ratio's) produce such different habitat rankings and how they relate to each other.   * If I can obtain log-ratio differences within the compana function.   * How I can view the data tables that compana has used. (I know how to view the tables that I entered myself, but `compana` replaced some zero-availability values by the weighted mean lambda, and I would like to see the new availabity table) ...Or any other advice on how to proceed.


title: MCAR assumption is plausible should I do MI?

text: My data consists of measurements on patients with cancer and the variables are some indicators regarding the cancer as well as the stage and the grade of the cancer and some personal info about the patients. According to the doctors the missing values occured either because the part of the tissue they took from the patients was destroyed or because the doctor who did the biopsy didn't take a tissue sample big enough to estimate the stage of the cancer. Taking all these into account I have no reason to believe that the variables are MNAR and so two options remain, MAR or MCAR. To my question now, let's say the MCAR is true (I read in another question that I can check this by doing a t-test, if someone knows another way please let me know) an easy way to deal with the missing values is to follow the complete-case method since the results will be unbiased, but my problem is that my data isn't big enough and if I delete cases I am afraid that the results probably won't be efficient. That's why I plan to do the MI but I'd like to know others opinion on that.


title: How to select cut-point for making classifications table for logistic regression?

text: How can I select cut-points to convert predicted probabilities to predicted responses in order to make a classification table for logistic regression? Should I take different cut-points like .5, .6, .7, etc.? (If I take different values the prediction error rate varies for them.) How can I make a generalization when I am taking different cut-points and having different prediction error rates?


title: Fitting multiple linear regression in R: autocorrelated residuals

text: I'm trying to estimate a multiple linear regression in R with an equation like this:               regr <- lm(rate ~ constant + askings + questions + 0)      askings and questions are quarterly data time-series, constructed with `askings <- ts(...)`. The problem now is that I got autocorrelated residuals. I know that it is possible to fit the regression using the gls function, but I don't know how to identify the correct AR or ARMA error structure which I have to implement in the gls function. I would try to estimate again now with,               gls(rate ~ constant + askings + questions + 0, correlation=corARMA(p=?,q=?))      but I'm unfortunately neither an R expert nor an statistical expert in general to identify p and q. I would be pleased If someone could give me a useful hint. Thank you very much in advance! Jo


title: How can I test the difference between a population proportion and sample proportion?

text: A report says that 82% of British Columbians over the age of $25$ are high school graduates. A survey of randomly selected residents of a certain city included $1290$ who were over the age of $25$, and $1012$ of them were high school graduates. Is the city's result of $1012$ unusually high, low, or neither? How should I approach this question? Any help and hints would be appreciated.


title: How can I do a correlation between Likert scale and an ordinal categorical measure?

text: I have a 5 pt Likert scale and a measure of time on the internet e.g. 1= less than 10 mins, 2=10-30 etc as well as number of SNS contacts (1=10-30, etc). I would like to estimate their correlation but I don't think I could do that directly as such. I computed a mean score for the Likert responses but what about the other two measures??


title: Max and min value of coefficient of determination

text: I have two questions: 1) What is the maximum and minimum possible values of $R^{2}$ for fitting a straight line? and 2) List all possible cases such that $R^{2}$ achieve max/min . I know $0\leq R^{2}\leq 1$. So max = 1 and min = 0. But how can i find all possible cases that makes $R^{2}$ achieve max or min? thanks


title: Strong correlations masking weak correlations?

text: I want to calculate correlations among thousands of genes using expression data across many different samples. For the non-biologists, I will try to explain. Basically, I have genes $i=1..M$ and samples $j=1..N$ and a $(M \times N)$ matrix $X$, with entries $X_{ij}$ corresponding to the expression levels of gene $i$ in sample $j$. Therefore, the correlations I am computing are between vectors $X_{k} = (X{_k}{_1}, X{_k}{_2},.., X{_k}{_N})$ and $X_l = (X{_l}{_1}, X{_l}{_2},..,X{_l}{_N})$. My problem is that sometimes strong correlations are masking weaker correlations that are still biologically relevant. One reason this may happen is that genes are organized into large functional modules that are themselves highly correlated (or anti-correlated). For example, maybe genes $X_5$ and $X{_1}{_0}$ belong to mutually exclusive modules and are therefore strongly anticorrelated. This anticorrelation masks weaker within-module correlations that may not be as significant when looking across all samples. The problem is I do not always know what functional modules each gene belongs to. I am looking for a general way to detect weaker correlations, perhaps by "subtracting" out the strong correlations in my data? Ideally, I would like to find many levels of correlation. This sounds like a simple statistical problem, but I'm not sure where to begin looking. I've looked at multi-level analysis, but I don't want to specify "levels" in my data. I am thinking more along the lines of "eigencorrelations" where each group of correlations describes different types of variation among the data. If anyone can point me in the right direction, it would be really helpful. Thanks!


title: explaining an extremely large coefficient in a rare events logistic regression

text: I am running a rare events logistic regression on a binary dependent variable. I have 538 observations and only 10 events (so 528 values of 0 and 10 of 1), which is why I chose to use a rare events logistic regression. When I run the regression, one of the independent variables in the model has a huge coefficient (around 25,000,000) and is found to be significant. The range on the independent variable is 0 to 1. Is this a problem? Could anyone explain why this is happening? When I run the same model with just a logistic regression this variable is insignificant. I'm not sure what is happening. Any advice would be appreciated.