title: What am I measuring when I apply a graded response model to the "Hunting of the Snark" dataset?

text: In another question, I asked about the statistical validity of StackExchange's "Hunting of the Snark" dataset, and whether or not we could draw any conclusions from its results. I measured a few reliability coefficients to better understand them; now I'm trying to create an item-response model, per Andy W.'s suggestion. First, a warning: I am not a qualified statistician in any discipline. It's been years since I ran a statistical analysis for any purpose, and at the time I wasn't terribly good at it. All the same, this is a lot of fun. I'm happy to run any analysis I can reasonably complete in R. The dataset is available here. I truncated all but the last 20 columns. I also cleaned the data slightly. In some cases, "Unfriendly" was misspelled as "UIfriendly". In row 4352, two respondents, 19 and 20, had '0' in place of a response. Given the way the rest of the data set is organized, I assumed this is null data, thus I nulled out those two responses. After reading this example of ordinal data modeling (specifically section 3.2) and this presentation about the 'ltm' package in R, I concluded my results were polytomous and I needed a graded response model. I first collected descriptive statistics of the data set. They suggested nothing surprising: "Neutral/unclear" was the most common response; no comment was universally friendly; some were universally unfriendly; respondents varied a fair bit in their frequency of neutral responses (see my other question for some measures of respondent agreement). My sources then suggested calculating the data's nonparametric correlation coefficients, which I have done here (you'll want to toggle text wrapping), using Kendall's tau. If I'm reading these results correctly, I believe they suggest some association between respondents. Finally, I ran the graded response models. One of my sources suggested fitting both a "constrained" and "unconstrained" model and comparing the two with an ANOVA test. Full results are linked to their respective models, below. I've pasted the summary lines here, as well as the results of the ANOVA: Constrained GRM:               Model Summary:         log.Lik      AIC    BIC         -96657.64 193401.3 193696      Unconstrained GRM:               Model Summary:         log.Lik      AIC    BIC         -96141.05 192406.1 192831      ANOVA:               Likelihood Ratio Table               AIC    BIC   log.Lik     LRT df p.value     fit1 193401.3 193696 -96657.64                        fit2 192406.1 192831 -96141.05 1033.18 19  <0.001      If I'm reading that ANOVA LRT correctly, it looks like the unconstrained model is a much better fit, so let's graph it. Everyone loves a graph. Out of the many I could post, I selected four: The combined Item Information Curves for the 20 respondents, and the Item Response Category Characteristic Curves for "Friendly" (category 1), "Neutral/unclear" (category 2), and "Unfriendly" (category 3). I'm happy to post any others that might be useful. They are, respectively: **Plot 1:** ![Combined IIC](http://i.stack.imgur.com/4vM1A.png) **Plot 2:** ![Category 1 -- "Friendly"](http://i.stack.imgur.com/Y2KyN.png) **Plot 3:** ![Category 2 -- "Neutral/unclear"](http://i.stack.imgur.com/f3SMA.png) **Plot 4:** ![Category 3 -- "Unfriendly"](http://i.stack.imgur.com/kKtVX.png) I lack the training to interpret these results, so I'm hoping someone out there can. What did I measure and what, if any, conclusions can I draw?


title: A summary of different cases of interaction effects, how do you interpret?

text: I know there are lots of questions/answers about dealing with interaction effects in regressions on this site, but I think there is a need to summarize a bit. I don't know if there is a standard way to deal with it, or if it makes sense. I just want to make a list of occasions when people may see interaction effects, and discuss about how to deal with them in each condition. If we are doing regressions with one covariate/independent variable/predictor or whatever you may call it, things are always easy. You just need to built a model like                              Y ~ B0 + X + error      and there is no need to bother interactions. If we have multiple covariates, say 2, things become complex. What I usually do is to build a model without interaction effects first:                              Y ~ B0 + X1 + X2 + error       Let's call it RWRTI (regression without interaction terms). Then I decide what to do next based on the results. Perhaps I build another model with interaction terms:                              Y ~ B0 + X1 + X2 + X1*X2 + error      Let's call it RWI (regression with interaction terms). If I firstly run RWRTI, I will list some cases with questions: _**(1) Both X1 and X2 are not shown to be significant in RWRTI._** For me, I simply report the non-significance of both covariates, and finish. I don't do RWI. I don't think anything would change with RWI. Do you think it might change? _**(2) Both X1 and X2 are shown to be siginificant in RWRTI._** For me, a natural reaction is to use RWI to check interaction terms. But this is sometimes problematic for me, e.g. After seeing the result from RWI, I might have 5 conditions :                1. X1 is significant, X2 is not, X1*X2 is significant.              This perhaps indicates that the significance of X2 is due to the interaction.          So I report X1 as significant, and for X2 I report the interaction effect.           2. Both X1 and X2 are significant, X1*X2 is not significant.              This is simple, since I just need to report the main effects of both.           3. All X1, X2 and X1*X2 are significant.              I need to report 2 main effects and an interaction           4. Both X1 and X2 are not significant, but the interaction term X1*X2 is significant               In this case, it means both significance of X1 and X2 are due to the interaction.           So I only report the interaction.           5. Nothing is significant.              I have never seen this happen.      _**(3) X1 is shown to be significant, but X2 is not with RWRTI._** For this case, I am a bit confusing. Since only one covariate is significant, should I still try RWI to see the interaction? Is it possible that the significance of X1 is due to its interaction with the other non-significant covariate X2? Before I thought I might be able to remove X2 from the model, since it is non-significant, but someone suggests not to remove it, because we cannot build models solely based on significance. Anyway, if I do not remove X2 and do a RWI, I will again facing the following conditions:                1. X1 is significant, X2 is not, X1*X2 is significant.              This perhaps indicates that the significance of X1 is partialy due to the interaction.          So I report X1 as significant, and I also report the interaction effect between the 2.           2. Only the interaction term X1*X2 is significant               In this case, it means the significance of X1 is from the interaction.           So I report the interaction effect only.           3. Nothing is significant.               This is something weried, but I really have this problem recently (for a mixed linear model).           The RWI model cancels out everything, no main effects, no interactions.           I have difficulty interpreting the results, because nothing is significant.           4. Both X1 and X2 are significant, X1*X2 is not significant.              I have never seen this happen.           5. All X1, X2 and X1*X2 are significant.              I have never seen this happen.      OK, I've listed all possible cases one might see when dealing with multiple regression. How do you deal with each case? Do you have any suggestions?


title: How to choose the right number of parameters in Logistic Regression?

text: I am studying Andrew Ng's Machine Learning lecture notes. I understand either we can manually choose the number of parameters, or we can use regularization to make it correctly fit. ![enter image description here](http://i.stack.imgur.com/XadCK.jpg) I was wondering are there any basic rules for choosing the right number of parameters? Can anyone give any explanation please?


title: Multiple hypothesis ANOVA

text: I'm not sure that I've titled this question correctly, but here is my query. Suppose you are given a set of measurements and the uncertainty (variance) associated with each. The task is to statistically figure out how many different objects were likely measured and finally, to combine measurements into a single estimate for each. The second part is easy enough - an uncertainty-weighted mean would do it - but I am having difficulty understanding how to sort out how many objects were measured. If there were just two objects, ANOVA would work. But what if there are an unknown number of objects? As an aside, I'm aware of a Baysian technique in which one considers each measurement in turn, building a set of hypotheses for each measurement that doesn't fall within the confidence interval of an existing hypothesis and combining it into the hypothesis when they do. But I think this method is dependent on the order the measurements are considered and therefore imparts a kind of time dependence on measurements that have none. I feel like this is something that's commonly done and I should know how to do, but I'm stumped so any thoughts you all might have would be much appreciated. Thanks, Val


title: Simulate data by using existing dataset

text: I have a complete dataset with input variables and response variables. I would like to perform a simulation where I give the input variables and generate randomly the response variables. Is there a way to do that without using parametric models (like in machine learning), or do I have to find a conditional joint distribution? I can already "predict" the value of one of the output variables by using random forest, but I can't do that for all the variables simultaneously.


title: Unbiased Estimator of Days Until Completion?

text: I'm trying to get an estimate of average number of days until some event occurs (the event is **guaranteed** to eventually occur). I have some sample where this event has already occured for most observations, so the number of days is known. However, this event has yet to occur for many other observations. I can't just omit those observations, since this will bias the average downwards. Can I get an unbiased estimate of the average? If so, how? Thank you!


title: Calibrating a multi-class boosted classifier

text: I have read Alexandru Niculescu-Mizil and Rich Caruana's paper "Obtaining Calibrated Probabilities from Boosting" and the discussion in this thread. However, I am still having trouble understanding and implementing **logistic** or **Platt's scaling** to calibrate the output of my multi-class boosting classifier (gentle-boost with decision stumps). I am somewhat familiar with generalized linear models, and I think I understand how the logistic and Platt's calibration methods work in the binary case, but am not sure I know how to extend the method described in the paper to the multi-class case. The classifier I am using outputs the following:   * $f_{ij}$ = Number of votes that the classifier casts for class $j$ for the sample $i$ that is being classified   * $y_i$ = Estimated class  At this point I have the following questions: **Q1:** Do I need to use a multinomial logit to estimate probabilities? or can I still do this with logistic regression (e.g. in a _1-vs-all_ fashion)? **Q2:** How should I define the intermediate target variables (e.g. as in Platt's scaling) for the multi-class case? **Q3:** I understand this might be a lot to ask, but would anybody be willing to sketch out the pseudo-code for this problem? (on a more practical level, I am interested in a solution in Matlab).


title: Optimal way to creating new features from training set (in R)

text: My training dataset has 1500 features - all numeric. (But only about 200 data points). I want to create additional features and then use the exhaustive list for feature selection. For creating additional features, I considered generating all 2-way interactions between the features. I ran the following code in R - but it is taking a lot of time. For 1500 features, selecting all 2-way products will result in about 1.2 million features. After about 2 hours of running, the result had only about 300k features generated. Is there a more efficient way to generate this? Could this be done faster? Here is the relevant code:               for(i in 1:ncol(train)){         for(j in i:ncol(train)){             two_way_features <- cbind(two_way_features,train[,i]*train[,j])     }}      EDIT: Apologies - I didn't add in the context. We are running an online campaign - sent to 200 companies. 6 measures were observed post the campaign. I have about 1500 features (including both about the campaign and about the companies). The aim is to learn from this. (And unfortunately, in B2B space, there's only so much we can do to get data. While I can try to get more features, can't expand to more than 200 customers/companies). Learning from this data, we would like to predict the 6 measures for another set of 100 companies (on whom we plan to run the campaign post thanksgiving). The entire dataset is numeric(both the measured target variables and the feature vectors). All of them are non-negative. No categorical variables too. For now, glmnet is what I am doing. After playing with glmnet for a while with the given feature vectors, the cross-validation score hasn't been improving. I am trying to investigate if combination of feature vectors provide a better CV score. Given that there's really limited data points, I haven't explored outside glmnet. (strictly LASSO).


title: Gibbs sampling from the full conditional distribution using R

text: I have a likelihood function of all the data $y$ $$L(\tau ,\theta |y)\propto \theta ^{\sum \delta _{i}^{C}+\sum \delta _{i}^{H}}\tau ^{\sum \delta _{i}^{H}}e^{-\theta \sum x_{i}^{C}-\tau\theta\sum x_{i}^{H}}$$ where $\sum \delta _{i}^{C}, \sum \delta_{i}^{H},\sum x_{i}^{C}$ and $\sum x_{i}^{H}$ are constant, and a conjugate prior $$f(\theta, \tau) \propto \theta^{a} \tau^{b}e^{-\theta c-\tau \theta d}$$ where the hyperparameters are $(a,b,c,d)=(3,1,60,120)$. I want to simulate samples from this distribution, but I have trouble in find the posterior distribution and full conditional distribution for the parameters $\tau$ and $ \theta$. From the prior distribution, can I infer that $f(\theta)$ follows $Gamma(a+1, c)$ and $f(\tau |\theta)$ follows $Gamma(b+1, \theta d)$? Hopefully someone can give me some references. Thanks very much.


title: $E(\frac{1}{1+x^2})$ under a Gaussian

text: This question is leading on from the following question. http://math.stackexchange.com/questions/360275/e1-1x2-under-a-normal- distribution Basically what is the $E\left(\frac{1}{1+x^2}\right)$ under a general Gaussian $\mathcal{N}(\mu,\sigma^2)$. I tried rewriting $\frac{1}{1+x^2}$ as a scalar mixture of Gaussians ($\propto \int\mathcal{N}(x|0,\tau^{-1})Ga(\tau|1/2,1/2)d\tau$). This also came to a halt, unless you folks have got a trick under your belt. If this integral isn't analytical any sensible bounds?


title: What do I need in order to draw conclusions from this data?

text: I have three techniques, called A, B and C. Each can be used independently when trying to perform four related tasks (Tasks 1, 2, 3 and 4). I have run lots of tests, trying all combinations of each technique being on or off. My results look something like this. Each number represents how many times the task was completed successfully when attempted, using the given combination of techniques, 100000 times. So, the higher the number, the better. $$ \begin{array}{l|r|r|r|r|r|r|r|r|r|} \mbox{Technique $A$} & - & - & - & - & X & X & X & X \\ \mbox{Technique $B$} & - & - & X & X & - & - & X & X \\ \mbox{Technique $C$} & - & X & - & X & - & X & - & X \\ \hline \mbox{Task $1$} & 433 & 277 & 911 & 492 & 686 & 4211 & 3775 & {\bf 9732}\\ \mbox{Task $2$} & 149 & 1063 & 5562 & {\bf 6035} & 3 & 58 & 1391 & 1708\\ \mbox{Task $3$} & 220 & 1278 & 7014 & {\bf 7018} & 10 & 97 & 2083 & 4452\\ \mbox{Task $4$} & 218 & 1255 & 6142 & {\bf 8656} & 1 & 73 & 1087 & 2056\\ \end{array} $$ I've highlighted the largest number in each row. Looking at the numbers, it seems that $B+C$ is good for Tasks 2, 3 and 4, and that adding $A$ as well is best for Task 1. But I want to say a bit more. I'd like to be quantitative if I can. My question is: can I deduce _anything_ quantitative from this data? Or do I really need some measure of the variance of the observations? That is, I suspect the numbers might be different if I ran all the tests again. * * * **Edit.** I originally presented my question pretending there were 3 techniques, just in order to save space. There are in fact 4. I'm mentioning this now because it might affect how the results are analysed. In particular: note now that the `intercept' is zero, i.e. there is no success at all when none of the techniques are used. $$ \begin{array}{l|r|r|r|r|r|r|r|r|r|} \mbox{Technique $Z$} & - & - & - & - & - & - & - & - & X & X & X & X & X & X & X & X \\ \mbox{Technique $A$} & - & - & - & - & X & X & X & X & - & - & - & - & X & X & X & X \\ \mbox{Technique $B$} & - & - & X & X & - & - & X & X & - & - & X & X & - & - & X & X \\ \mbox{Technique $C$} & - & X & - & X & - & X & - & X & - & X & - & X & - & X & - & X \\ \hline \mbox{Task $1$} & 0 & 0 & 0 & 0 & 0 & 2700 & 1 & 10113 & 433 & 277 & 911 & 492 & 686 & 4211 & 3775 & {\bf 9732}\\ \mbox{Task $2$} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 149 & 1063 & 5562 & {\bf 6035} & 3 & 58 & 1391 & 1708\\ \mbox{Task $3$} & 0 & 0 & 0 & 0 & 0 & 664 & 0 & 3043 & 220 & 1278 & 7014 & {\bf 7018} & 10 & 97 & 2083 & 4452\\ \mbox{Task $4$} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 218 & 1255 & 6142 & {\bf 8656} & 1 & 73 & 1087 & 2056\\ \end{array} $$ Here are those numbers in by-column form (for pasting into `R`):               0,    0,    0,    0,     0,    0,    0,    0,     0,    0,    0,    0,     0,    0,    0,    0,     0,    0,    0,    0,     2700, 0,    664,  0,     1,    0,    0,    0,     10113,0,    3043, 0,     433,  149,  220,  218,     277,  1063, 1278, 1255,     911,  5562, 7014, 6142,     492,  6034, 7018, 8656,     686,  3,    10,   1,     4211, 48,   97,   73,     3775, 1391, 2083, 1087,     9732, 1708, 4452, 2056


title: Diagnostic for VAR model. non normal

text: I have some problem about my model. my model is based on VAR. (vector auto-.) well, I've tested ARCH test, BG test(autocorrelation p) and jarque.bera.test. Model is stable. Also I got good result for ARCH and BG (both fail to reject) but I failed jarque.bera.test. What does it mean? I know it didn't satisfies white noise assumption but I was wondering if this suggest other model to use. Or is this because of omitted variable problem?


title: Transformation of the datasets with negative numbers for exponential graph?

text: I have a simple question for data transformation for fitting my dataset to a negative exponential graph. There are negative values on my dataset which hinder fitting a negative exponential curve. How can I successfully transform my data in statistical way? My data is about Mantel correlation coefficients with distance. Here is my data.               Distance   Coefficients         1       0.7241232         2       0.19643728         3      -0.06509062         4      -0.16492022         5      -0.39367865         6      -0.50758682      What is the best way to transform this data?


title: Book on interpreting plots in regression analysis

text: I am taking a regression course that is mainly theoretical. But we are given scatter plots and residual plots, and we are supposed to be able to describe the relationships and interpret them. I find interpretation part difficult sometimes. Is there a good book that has exercises using SPSS output that I can use to help me interpret regression plots?


title: SVMs and separating margins

text: A couple of points that I need cleared up, thanks ahead for any help.   * In a D dimensional SVM what is the dimension of the maximal separating plane? Intuitively (and following VC logic) it will be D-1, is that correct....   * In the derivation of the SVM, what roles do the Lagrangian multiplier play? What do their values indicate?


title: Exponential family in testing and estimation

text: In the Annals of Statistics paper "Defining the curvature of a statistical problem(with applications to second order efficiency)" by Bradley Efron, he claims the following two statements in the first paragraph.   1. The locally most powerful test of $\theta=\theta_0$ versus $\theta>\theta_0$ is uniformly most powerful in an exponential family.   2. The MLE for $\theta$ is a sufficient statistic in an exponential family and achieves the Cramer-Rao lower bound if we have chosen the right function of $\theta$ to estimate. Can someone elaborate on these two statements as I am a Mathematician who has only a little knowledge of Statistics?


title: Run time analysis of the clustering algorithm (k-means)

text: I was reading some notes on ML and clustering and it claimed that the run time of clustering was O(kn) where k is the number of clusters and n is the number of points. I was wondering why this was true and if someone had an analysis for it. This is what I have thought so far though: Recall the clustering algorithm. 1) Initialize centroid $z^{(i)}, ..., z^{(k)}$ 2) Repeat the following two steps until there is no further change in cost: a) For each fixed centroid $z^{(i)}$, choose a optimal clustering such that point $x^{(i)}$ is closest to the current centroid. b) For a fix cluster, re-compute the current centroids: $z^{(j)} = \frac{\sum_{i \in C_{j}} x^{(i)}}{|C_j|}$ We know that clustering does not run forever because of the convergence proof. So it doesn't run forever. Now, step a) takes O(kn) worst case because for each centroid $z^{(i)}$, we have to compute the Euclidean distance to each data point $x^{i}$ and choose the min one and assign it to that cluster which in the worst case we do O(n) to compute min for each centroid k. So O(kn) (I know we can use heaps or whatever, thats not what interests me). Now on the next step we compute the new means which worst case takes O(k+n) which is less than O(kn) so step b) doesn't really contribute anything to the runtime. However, Basically...how do we know that it won't keep doing this for a lot more times? How is it that the notes imply that it converges in O(1). That part was not clear to me. If its clear to anyone else, feel free enlighten me!


title: Distribution for time between arrivals at the same location by the same entity

text: Given a sample of inter arrival times for an entity at the same location, is there a distribution that can fit such a sample?. I was thinking about poisson, but was concerned about the independence, since the next arrival may be dependent on the previous arrival time. Or is using a time series approach better to capture trend and seasonality?.


title: Non-Parametric Regression with Ordinal Independent Variables

text: I'm trying to find a model to predict a test score (0 to 100) from many independent ordinal variables (all are on a 1 to 10 scale). Since there are SO many ordinal variables, using binary dummy variables would probably be a mess. The test scores are skewed, (as well as many of the ordinal variables) so I'd like to use a non-parametric approach. Another issue/thing to keep in mind is that a score of 10 on the predictor variables does not necessarily predict a higher test score. For example: a score of 2 in one predictor may predict a high test score, while a score of 8 on another may predict a high test score. Not sure what model to use and how to deal with all of the ordinal variables. I'm working in R, so please put any coding advice in R. Thank you!


title: Correlation analysis and correcting p-values for multiple testing II

text: I am currently facing problems dealing with lots of individual simple correlations (28 individual DTI connections strengths, individual behavioural scores on 31 time points). That would be 28x31 individual correlations. Referring to a previous discussion on that, I just want to make sure that I understood correctly and am looking forward for your short answers whether this way would be fine:   1. Individual correlations 28x31 are computed, giving R and p values (original data).   2. Subjects are shuffled on 1 side (e.g. for all 28 connections, pseudosample), all correlations are computed again, giving a new 28x31x5000 matrix with R and p-val. 2a. This is done 5000x.   3. The number R perm > R obs (N) is collected for each correlation (hence max. 5000).   4. Emperical p val is computed as N/5000 for each correlation.   5. A common threshold of p-val of e.g. 0.05 is applied to find significant correlations based on all individual 28x31 p-vals. My additional questions would be: A. What about negative R? What about the case when abs(R) would be larger in perm. than in original?   B. How to adress R perm -0.05 and R original -0.6, hence stronger negative correlation? At this stage, I only count the # of Rorig > Rperm, should I also count abs(Rorig> abs (Rperm)? I am still a bit confused and unsure whether I would report correct statistics.


title: Experimental design for testing changes in categorical outcomes

text: The setup: I am testing, if different instructions on playing a game result in different outcome. The game has fixed amount of N different possible outcomes. (I most likely set N to something around 2-5). The outcomes are measured on categorical scale (each possible outcome being one category). Lets say I have two different instructions: A and B. I want to plan an experimental design to test hypothesis: For each person, the outcome of the game is same for instructions A and B. Alternative hypothesis is: For each person, the outcome of the game depend on the instruction. The amount of test subjects could be around 50-200. (I prefer less for economic reasons) Do you have idea of what kind of experimental setup and statistical tests I can use?


title: Modeling a spatial trend by regression with the $(x,y)$ coordinates as predictors

text: I plan to include coordinates as covariates in the regression equation in order to adjust for the spatial trend that exists in the data. After that, I want to test residuals on spatial autocorrelation in random variation. I have several questions:   1. Should I perform linear regression in which only independent variables are $x$ and $y$ coordinates and then tests residuals on spatial autocorrelation, or should I rather include not only coordinates as covariates but also other variables and then test residuals.   2. If I expect to have quadratic trend, and then include not only $x,y$, but also , $xy$, $x^2$ and $y^2$, but then some of them ($xy$ and $y^2$) have the $p$-value higher than the threshold --should I exclude those variables with higher $p$-value as being nonsignificant? How should I then interpret the trend, it is certainly not quadratic anymore?   3. I guess I should treat $x$ and $y$ coordinates as any other covariates, and test them on having linear relationship with dependent variable by constructing partial residual plots ... but then once I transform them (if they show they need transformation), that will not be that kind of trend any more (especially if I include $xy$, $x^2$ and $y^2$ for quadratic trend). It may show that $x^2$, for example, needs transformation, while $x$ does not or so? How should I react in these situations? Thank you.


title: Under what conditions can a PLS regression model be expressed by single linear equation?

text: I am confused by two, yet inconsistent for me, facts: Since the PLS regression is expressed by matrices of scores and loadings as $$X=TP^T+E\\Y=UQ^T+F$$ how it can be translated into linear equation like $Y=a+b_1X_1+b_2X_2+...+b_nX_n$? (I have found this in several papers).


title: Which glm algorithm to use when predictors are numerical as well as categorical?

text: I just need a direction on which regression algorithm (preferably glm or similar) algorithm to use when the predictor variables are a mix of numerical and categorical variables. The output is numerical for the time being but in future, I need to extend this for categorial output also. My input columns are in the the format given below:               preNum1 predNum2 ......predCat1 predCat2 ....ResponseVar        Edit: I am trying to predict the amount of 'Bilirubin' (here, the response variable), depending on patients' data. We have patients's lab results such as glucose, blood pressure, etc etc as numerical predictor values and the patients' disease history (diabetes,and other diseases) as categorical predictors. I converted all the categorical variables (they are diseases such as diabetes, high BP, etc) into 0 and 1 (0 representing that it is not present, 1 representing it to be present.) After this I applied GLM on this. I am wondering if my approach is correct or not?


title: Can I categorize the factor scores to use them as predictors of an ordinal logistic regression?

text: I was wondering if I can categorize the factor saved scores by taking their quartiles (or some other measures, I am not sure what should I use!) as cut points and use them as predictors in an ordinal logistic regression. The reason I am thinking of doing so is that, the number of empty cells increases to a lot extent if I use the factor saved scores as they are (continuous in nature). I have no particular reason behind choosing quartiles as the cut points. The factor scores were obtained from an optimal scaling of the variables (some variable were continuous, some were ordinal, some were nominal) and then factor analysis of the optimally scaled variables using varimax rotation. I am just thinking if I can categorize the factor scores by some means and try to name them as, probably, "lowest scores", "lower scores", "higher scores" and "highest scores" etc. and use them in an ordinal logistic regression as independent variables so that I can comment on the different levels of scores compared to the reference category ("lowest scores") of that factor. I am not sure what should be used as the cut points. Even I am not sure if this is a right procedure. Can anyone please suggest? Thanks in advance for your kind theoretical support.


title: How to compute the Kullback-Leibler divergence when the PMF contains 0s?

text: I have the following timeseries ![enter image description here](http://i.stack.imgur.com/Q6Qld.png) obtained using the data posted below. For a sliding window size of 10, I am trying to compute the KL-divergence between the PMF of values within the current sliding window and the PMF of the history with the final goal of plotting the value of KL-divergence across time so that I can compare two time series. As of now, there is a conceptual problem I am facing (which I'll explain using Python):               In [228]: samples = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1]          # In reality this 10 should be 20 because that is the max value I have seen in the timeseries     In [229]: bins = scipy.linspace(0, 10, 21)      In [230]: bins     Out[230]:     array([  0. ,   0.5,   1. ,   1.5,   2. ,   2.5,   3. ,   3.5,   4. ,              4.5,   5. ,   5.5,   6. ,   6.5,   7. ,   7.5,   8. ,   8.5,              9. ,   9.5,  10. ])     In [231]: scipy.histogram(samples, bins=bins, density=True)     Out[231]:     (array([ 1.63636364,  0.        ,  0.36363636,  0.        ,  0.        ,             0.        ,  0.        ,  0.        ,  0.        ,  0.        ,             0.        ,  0.        ,  0.        ,  0.        ,  0.        ,             0.        ,  0.        ,  0.        ,  0.        ,  0.        ]),      array([  0. ,   0.5,   1. ,   1.5,   2. ,   2.5,   3. ,   3.5,   4. ,              4.5,   5. ,   5.5,   6. ,   6.5,   7. ,   7.5,   8. ,   8.5,              9. ,   9.5,  10. ]))      The problem is that the resulting `PMF` contains `0`s so that I cannot really multiple two `PMF`s to get the KL-divergence. Can someone tell me how to mitigate this problem? ### Data               0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 4 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 2 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0


title: small n: some variables are normally distributed some are not

text: We performed a biomechanical study, comparing an intact state to 5 different deficiency states on the same knee cadaver specimen. This was conducted on a very small sample size (n=8). My problem is, that some of the 5 deficiency variables are not normally distributed. My original plan was to compare each of the deficiency variable with the intact variable via a paired t-test. As some of the variables are not normally distributed my questions are:   1. Do I have to use a non-parametric test (Wilcoxon) for all variables or only for those, which are not normally distributed.   2. Can I still use paired t-test, when both variables are normally distributed   3. When I want to apply Bonferroni Correction, do I have to divide the Alpha by the number of total tests performed(parametric and non parametric) or only by the number of tests performed from the same category.


title: Leverage - influential points

text: Do we look at the absolute value of the leverage or the relative value? For instance, based on the chart below, the largest leverage is about 0.023, it is big compared to other data points, but I'm not sure if there's something like threshold in VIF that indicating a high leverage? ![Leverage Plot](http://i.stack.imgur.com/ISS0f.jpg)


title: Two fair coins, why increasing trials in prop.test not reduce false positive?

text: I am simulating two fair coins and running a `prop.test`, how come increasing number of trials doesn't decrease false positive? I'm suspecting it's in some underlying assumptions in the R functions that I'm using which I missed. This binomial coin flip simulator takes argument N, for number of trials, and returns the p-value of an two-sided `prop.test`.               test.flips <- function(N, A.prop=0.5, B.prop=0.5) {       heads.A <- rbinom(1, N, A.prop)       heads.B <- rbinom(1, N, B.prop)       test <- prop.test(c(heads.A, heads.B), n=c(N, N), alternative="two.sided")       return(test$p.value)     }      I do 100 times of 1000 trials each to get this set of p-values, i.e. `p.values <\- replicate(100, test.flips(1000))` ![enter image description here](http://i.stack.imgur.com/oga8O.png) The red horizontal line denotes a 0.05 p-value. Notice that we have a few false positives. Thus I do a power calculation, `power <\- power.prop.test(p1=0.5, p2=0.501, power=0.90, alternative="two.sided")` to find that `N = 5253704` if `power=0.90`. So I do this again, `p.values <\- replicate(100, test.flips(round(power$n)))` with N set to that high number of trials. But the number of false positives didn't improve as shown here: ![enter image description here](http://i.stack.imgur.com/VfHeK.png) If all else remains equal (same sig. level), how come I'm not seeing what's gained by increasing sample size? What's wrong and how do I fix this please? Edit:   * state clearer question in end   * As requested, here's `prop.test()$statistic` versus sample size. ![enter image description here](http://i.stack.imgur.com/GMBqH.png)


title: Probability of Default

text: I'm doing a project to predict probability of delinquent for individual loans. Seems the model I fit is not good and I want to improve the model. However, I'm confused by the results I got and don't know what to do next. Can anyone kindly give me some instructions on the project? The outcome is binary, where 1 stands for failing to pay in the next month and 0 stands for successfully making the payment. There are both categorical and continuous variables in the model. The data I have include loans originated within 1998 to 2007, with their original information, like loan size, zip code, credit score, interest rate, purpose code (primary residence or investing) etc. and dynamic information, like payment history, loan age, current fico score, current loan to value ratio, cumulative home price appreciation etc. The entire data table is too huge, so I did random sample selection stratified by origination year to get a sample dataset with 1 million loans. You can see my approach of treating data here -- although loan data is time series, I treat them as individual data points. Say loan A has a perfect payment history as of 1/1/2005 (that is the borrower made every payment since loan origination till now), I assume the probability of the borrower failing to make payment in 2/1/2005 solely depends on the variables I mentioned above. To estimate the coefficients of the variables, I ran a logistic regression. Since the model includes continuous variable, I think it is logistic regression with ungrouped data, so overdispersion wouldn't be a concern here. The following is what I get by analyzing the result in SAS. ![Delinquent Rate Actual VS Estimate](http://i.stack.imgur.com/qvT6Q.jpg) The chart above (actual weighted average delinquency rate vs predicted weighted average delinquency rate) looks fine, however, if we checked the pseudo R square, it’s only 0.0826. Log-likelihood ratio shows we should reject null hypothesis. Hosmer Lemeshow is significant, indicating the model doesn't fit well. ROC is 0.64, seems prediction power is quite poor. I also used `proc glm` with tolerance option to test for multicollinearity, the code is below.               proc glm data = ltm.smp_PA_fixed_amort_3;     class purpose_code prop_type_desc fmonth FICO_bucket year_range(...);     model DLQ = orig_amt purpose_code FICO_bucket months_in_DLQ     loan_age prop_type_desc cumu_HPA fmonth year_range(...)/tolerance;     run;      However, I'm not sure how to interpret VIF of categorical variables, since each level of a categorical variables has its own tolerance. Last, I also checked the specification error. Used the code below:               proc logistic data = temp.PA_fixed_amort_3_dlq_5;     model dlq = pred pred2;     output out = temp.speci_test;     ods output parameterestimates = temp.speci_test_coeff;     run;      Here `pred2 = pred**2`, `pred` is the linear predictor from the logistic regression. Both `pred` and `pred2` are significant, indicating misspecification. Summary: It looks good when we looked at the fitted chart; A projection forward using Monte Carlo simulation also shows OK results for now (we project from 2011 to 2020 and compare with the actual rate from 2011 to 2014). But I doubt the robustness of the model, given many statistics indicating bad fit.   * Can anyone give me some thoughts about where to go next?   * Is there any issue with my approach (like identifying multicollinearity, specification error) or my interpretation of the results?   * How to utilize the tolerance data to decide if a categorical variable is correlated with other variables / should be removed or transformed?


title: calculating standard error of coefficients for distance-based permutation model in R

text: I apologize that this seems to be a half function-specific, half stats question. I'm conducting a set of multivariate analyses using the adonis function from the vegan package in R - distance-based permutational manova. Predictor variables are mixed continuous and categorical. I'm using the same predictors to fit 3 models to 3 sets of response variables - each set being all binary, ordinal, or continuous (count) data. I would like to display results by plotting the coefficients with confidence intervals, as discussed here: tables2graphs.com/doku.php. (The approaches outlined at those links would need to be modified for multivariate models - that's not my concern.) I've searched a good chunk of several days on the following questions, without much luck: (1) Does it make sense to even consider coefficient SEs with a permutational distance-based method such as adonis? (2) (If so) How could I calculate coefficient SEs from the adonis output and/or the raw data? Looking at this formula for coefficient SE (found here: http://stattrek.com/regression/slope-confidence-interval.aspx) SE = sb1 = sqrt [ Σ(yi - ŷi)2 / (n - 2) ] / sqrt [ Σ(xi - x)2 ] "where yi is the value of the dependent variable for observation i, ŷi is estimated value of the dependent variable for observation i, xi is the observed value of the independent variable for observation i, x is the mean of the independent variable, and n is the number of observations." I run into confusion with the mixed nature of the data, and come up with the following questions: (1) Could I use the term-level MSE from the anova table for the first half of that equation? Otherwise (if it would need to be recalculated for each factor level) I'm stuck at adonis' lack of a predict function. (2) If so, can the second half of the equation be calculated using the coding in the model.matrix for the factors? Or is R doing something sneaky with the coding that I don't understand, rather than using it in as a set of values in a straightforward way? Any help or guidance is greatly appreciate. Thank you for your time and consideration! Warmly, Rafter


title: Non negative least squares with minimal colinearity

text: I am trying to fit a dataset using the standard NNLS (non-negative least squares) approach. Formally: $\min_x ||Ax-b||^2_2$ s.t. $x\ge0$ This is a quadratic program and can be solved optimally. The solution I find fits the data reasonably well and is relatively sparse (which is good for me) but has an undesirable property: I find that it may assign high weights to features (columns of $A$) that are highly correlated. I would like my solution to be such that the correlation between support features (features that have high weights) will be minimal. Note that in my setting all entries of $A$ and $b$ are non-negative, and I can normalize the columns of $A$ so that the norm of each column is 1. I tried approaching the problem by directly minimizing possible formalizations of the quantity I want. Note that $A^TA$ is the covariance matrix, so it could potentially be a good thing to minimize (perhaps not the diagonal, though). But minimizing $||Ax-b||^2_2+x^TA^TAx$ does not make any sense and minimizing $||Ax-b||^2_2+x^T(A^TA-I)x$ makes the optimization non-convex (because you get negative eigenvalues). With this approach, I can see that intuitively this is a bit like doing the opposite of ridge-regression/Tikhonov regularization. I also tried L1 regularization ($||Ax-b||^2_2+\lambda||x||_1$) just for the heck of it, but it doesn't solve this problem. Formally, I am looking to modify the objective function such that it penalizes solutions in which features that are "similar" to each other both get high weights. There are several ways to formalize this notion, for example one being that I would like the support columns to be as orthogonal as possible. Does anyone have ideas on how else to approach this?


title: Granger causality

text: How can you calculate the granger causality between x(t) and y(t) when you have only been given two equations and no data? I have been stuck on this for a while. Any help will be really appreciated.


title: Should PCA be performed before I do classification?

text: I have got a problem about doing a classification. I have got around 50 datasets. Each of them has 15 features. I am trying to use these features to classify the 50 datasets to either 'Good' or 'Bad'. The ground truth labels of the 50 datasets are available so that a classical training and validation can be done. As there are 15 features, the problem should be considered as a classification in high dimensions. My question is: Should we always perform PCA before we run any generic classification algorithms, such as LDA, KNN or SVM? I got someone's opinion that: "PCA chooses the directions in which the variables have the most spread, not the dimensions that have the most relative distances between clustered subclasses." But for my understanding, in order to do a better classification, we need to find features that have large differences between two groups. For example, we can calculate the mean and standard deviation of a feature for 'Good' and 'Bad' separately, and we can see if there are large differences. If so, we choose this feature. Also, we need to find features that have got least correlation in between. If two features have a large positive correlation, we can just choose to use one of them. PCA is somehow picking up the dimension reduced features for us, given 15 features, it will give 2 or 3 principal components that can be classified better. Am I right? Or I am on the wrong course?


title: MCA to Predict New cases?

text: Just wondering if we can use MCA (Multiple Correspondence Analysis) to predict new cases's classification?   1. Map new cases on to the Biplot to see which direction on the axes it is closer to?   2. How do we calculate the distances of the new cases? Is there a way to implement thess algorithms in R? Thanks!


title: Unsupervised feature learning from raw text as a previous step for clasification?

text: I have a corpus of 2500 opinions, is it posible to use scikit´s restricted boltzmann machine implementation to extract a feature vector as a previous step to a classification task?. What aproach do i need to follow in order to use a restricted boltzman machine for text classification?, do i need to label my data to classify?


title: How to prove that the manifold assumption is correct?

text: In machine learning, it is often assumed that a data set lies on a smooth low- dimensional manifold (the manifold assumption), but is there any way to prove that assuming certain conditions are satisfied, then the data set is indeed (approximately) generated from a low-dimensional smooth manifold? For example, given a data sequence $\{\mathbf{X}_1 \ldots \mathbf{X}_n\}$ where $\mathbf{X_i} \in \mathbb{R}^d$ (say the sequence of face images with different angles) and a corresponding label sequence $\{ y_1 \ldots y_n\}$ where $y_1 \preceq y_2 \ldots \preceq y_n$ (say the angles of the face sequence). Suppose when $X_i$ and $X_{i+1}$ are very close, their labels $y_i$ and $y_{i+1}$ are also very close, we can imagine that it is likely that $\{\mathbf{X}_1 \ldots \mathbf{X}_n\}$ lie on a low-dimensional manifold. Is this true? If so, how can we prove it? Or what conditions does the sequence need to satisfy so that the manifold assumption can be proven to be true? Thanks in advance!


title: Imputing values that depend on each other, such as percentiles or proportions

text: I have a data set with school level measures including test score percentiles. These percentiles are central to my analysis. For example, one measure is the test score of the 25th and 75th percentile student at that school. Likewise I have data on race that needs to sum to 100. I have missing data and want to use multiple imputation. These values depend on each other because, eg, a 25th pct score can't be higher than a 75th. Likewise, if 90 percent of a school is white, 25 percent can't be black. I do not have the individual level data at each school upon which those measures are based. I do however have individual test scores of students, just not enough at each school to estimate it. In MI procedures I've run so far I have a problem with some imputed cases not having realistic values. In some schools it imputes 25th percentile scores that are larger than 75th. Likewise, the percents of various races do not add up to 100. I could write some algorithms that adjust the imputed values to be realistic (eg shrinking or inflating percents proportionally until they add to 100). This seems ad hoc and hackish. One option is to use the individual level measures I have to fill in plausible values. For example, I have test scores on individuals so I know the typical 25tg pct score. I can also adjust these for different school types (eg public or private etc). This seems problematic in that there will be too little variation among imputed cases, implausibly lowering variance. It may also make sense to use hot deck imputation but I've never done this before. Anyone have tips? I'm stuck between a lot of imperfect options. I'm leaning toward using my individual level data to produce realistic values that account for some school differences. But I'm wondering if anyone has any ideas on how to responsibly use MI. I use stata with imputation by chained equations.


title: Estimating ability using IRT when the model parameters are known

text: I have 3PL model parameters (guessing, difficulty and discrimination item parameters). Is there any function with which I can estimate individual ability from item response data? I tried the function `factor.score` in the package ltm, but it seems to require the whole data from which the parameters were estimated, which I don't have.


library(ISLR)     lm.fit = lm(mpg[-1]~weight[-1], Auto)     predict(lm.fit, data.frame(weight = Auto$weight[1]))      Why does this code cause the following error:               Warning message:     'newdata' had 1 row but variables found have 0 rows      It doesn't cause an error for any other index but 1.


title: Coding for an ordered covariate

text: I am performing ordinal regression, I have 5 response categories and several predictors both continuous and categorical. I would like to add a predictor which is categorical but ordered (1, 2, 3, 4). I don't think it would be appropriate to apply the usual dummy coding for unordered categorical predictors, but when I searched for how to code this I did not find much information. In Steyerberg (2009) "linear coding" or "assuming linearity of the predictor effect" is mentioned, but without further details. Does it mean I just use my ordered values as they are, i.e. use them as a continuous variable?


title: Creating Bin in graph by R from a set of data

text: I have 1000 of data of mean temperature and daily rainfall of specific area. I need is to make a graph by using R software. I had already made the structure and summary of those given data. Now I have to make bins of 10 sets of temperature data and determine the 95th percentile rainfall of each bins. Next is to do the log of each bin's 95th percentile rainfall data in y-axis and the mean temperature of each bin in x-axis. How will i do it?


title: Omitted Variables because of collinearity

text: memb            float   %9.0g                 proportion of permanent employees who are                                                     members of the worker cooperative     bonus           float   %9.0g                 average distributed profits per worker in                                                     millions of lire          . reg lnQ lnL lnK memb own lnBonus     note: memb omitted because of collinearity     note: lnBonus omitted because of collinearity                Source |       SS       df       MS              Number of obs =       4     -------------+------------------------------           F(  3,     0) =       .            Model |    2.353812     3  .784603999           Prob > F      =       .         Residual |           0     0           .           R-squared     =  1.0000     -------------+------------------------------           Adj R-squared =       .            Total |    2.353812     3  .784603999           Root MSE      =       0          ------------------------------------------------------------------------------              lnQ |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------              lnL |   .4852686          .        .       .            .           .              lnK |   .4817975          .        .       .            .           .             memb |          0  (omitted)              own |   .0133421          .        .       .            .           .          lnBonus |          0  (omitted)            _cons |   2.260949          .        .       .            .           .      when I run a regression, memb and InBonud are omitted. I know it says because of collinearity....However, based on the description for both variables , I do not see how they relate to each other. Does anyone have an explanation why memb and InBonus are omitted?


title: Fisher's distribution and its approximation

text: I would like to ask how I can illustrate that Fisher's distribution when the null hypothesis is right goes to $F(k-m,n-k)$ when the degrees of freedom grow. How exactly is $F$ calculated in this case? And what is $k$? Should I use Wald test here? (It's assumed that errors have t distribution.)


title: Log-scaled chart for visualizing extreme range

text: I need a simple bar chart or the like. I am not a statistician by any means, but this chart is supposed to accurately represent FileIO in MB/s compared to the theoretical peak of a specific drive. The first value is always 30,000.   The second value is variable, but leans towards the extremely low 1–10 thousand-ish area. Can anybody recommend how to scale this the best way so that you can see changes in lower values and that they're really far from the peak?


title: rknn - Error: too many ties in knn

text: I am trying to run rknn with rknn package but it appears an error I am not able to understand. Here's the code and the error:               rknn(bank_training, bank_testing, y_n, k=2, r=100, mtry=4)          Error in knn(train = data[, fset], test = newdata[, fset], cl = y, k = k,  :      too many ties in knn      I tried to run it with different values of r from 10-500 but i get the same error and I don't understand what is meant with that. I ran it also as it is described in the manual (See below), but i get always the same error.               rknn(data, newdata, y, k = 1, r = 500, mtry = trunc(sqrt(ncol(data))),     Random.seed = NULL, seed = NULL, knn.algo="VR")      Does anyone have an idea why is that? Thanks in advance!


title: Metropolis random walk using multivariate normal

text: I need to implement a program that generates a sample from a really complicated distribution $f$ of 3 variables. I need to implement it using Metropolis-Hastings Algorithm's and its variation and I was thinking to use Metropolis Random Walk with a multivariate Normal distribution. But to use Random-Walk, the proposal must be symmetric. Is the multivariate normal always symmetric or is there any condition (about the matrix of covariances or anything else) that it has to satisfy so it is symmetric? Thanks in advance!


title: Testing hypothesis & $\alpha$

text: Hey I have a question that I found in a textbook for some practice before a test, but there aren't any solutions for it. I'm pretty sure it's related to testing a hypothesis, but I'm not sure. If anyone can point me into the right direction of solving this that would be great. Q1: Five soft drink bottling companies have agreed to implement a time management program in hopes of increasing productivity (measured in cases of soft drinks bottled per hour). The number of cases of soft drinks bottled per hour before and after implementation of the program are listed below.                       Company          #  1   2   3   4   5     Before 500 475 525 490 530     After  510 480 525 495 533      Test at $\alpha$ = .05 if the time management program is efficient in increasing the productivity.


title: Tarkkonen's rho - Reliability or Validity

text: Kindly advise what the value we obtain for Tarkkonen's rho tells us about - reliability or validity of the data being measured. Thank you


title: Interpreting a lots of effects

text: Does anybody know how to interpret a whole bunch of effects (main and interaction) in a clever way? Or does anybody have a good example where it's shown? To be more precisely: Assume that you have a lot of effects in your model (main and interaction effects) and you know that standard errrors and coefficients are biased (as they are in most empiric studies). But you want to present results, what would you do? So because of the bias p-values are wrong, so there is no valid way to say which coeficients are significant and which are not. So there has to be another way to present information. The idea to work with plots seems to be a good one. But I'm still curious how you guys handle such problems!


title: Modelling adwords clicks with Bayesian Stats

text: How can I use bayesian statistics for adwords clicks modelling? Is there any good source to learn about it?


title: Gibbs sampling with Log-Normal observations

text: I am writing a Gibbs sampler for data that is Log-Normal (LN) distributed, with unknown mean and variance. There is a wealth of information on inference for LN models when either the mean or variance (precision) are known, but I'm not finding much information on inferring both parameters. I have an idea of "reasonable" bounds on both parameters, but otherwise want to remain mostly uninformative. This question on priors for LNs seems related, but they simply mention that one should use Jeffrey's prior for the variance. I'm looking for more specific information: which priors and how would one sample from them. Also, would we have to resort to a Metropolis step or can we make some assumptions in the priors such that we can get posteriors in closed form?


title: Are significance tests only about sampling error or other forms of variation also?

text: When we study differences between a treatment and control group in an experimental setting, we want to test whether those differences are big enough that we can conclude they are not simply due to chance. _(Lacking official terminology for the some concepts below [especially #2], please either grant me a little creative license or suggest the appropriate technical terms)_ Differences due simply to chance, i.e., stochastic variation, could come from:   1. Sampling error, if the sample is smaller than the population.   2. Inherent error (or parallel universe error) -- There is some inherent unmeasurable randomness in the manifestation of an outcome in any given individual. There is a broader universe of possible outcomes than what we actually observe. Even if we measure the whole population, there is some fundamental level of variation where these observations could have been different, even if nothing meaningful or causal changed.   3. Measurement error -- Given imperfect questions, questioning, or responses, we are not measuring exactly the concept we are aiming for.    4. Temporal variation -- the moment at which we interview someone may lead to variation in answers even if they answer with perfect accuracy. e.g., "how much money did you spend on good X in the past 7 days?" This would vary depending on which day we interview the person. For each individual, we observe an outcome with some noise. $y_i = X_i \beta + \epsilon_i$ But the noise is due to a mix of the factors above. In my specific case, my team has surveyed 76% of our entire population of interest (due to 24% non-response rate). If we limit our scope conditions to the people who responded to the survey, then we have surveyed 100% and we have no sampling error. My colleague argues that we do not need significance tests. For example, if we observe that group A has $\bar{X_1} = 40$ and group B has $\bar{X_1} = 41$, we can infer that this is a meaningful difference. My perspective is that we should incorporate the variance of the measures because the magnitude of difference and variance of the data can still tell us about whether we are observing _systematic variation_ between groups as opposed to _stochastic variation_. Even though if we had measured our entire "real" population of interest, Due to stochastic variation from $\epsilon_i$ (or if $\epsilon_i$ is correlated with group assignment) we might observe differences between the two groups that we cannot interpret as "substantive," or "an effect of treatment" in an experimental setting. **Is my colleague correct that significance tests are exclusively about sampling error?** **How should a researcher test for stochastic variation other than sampling error?**


title: Fisher's & McNemar's tests

text: I'm conducting a retrospective study to assess the efficacy of a new test for predicting sudden cardiac death (SCD). I also want to compare the accuracy of my new predictor with that of a conventional predictor. I'm not a statistician, so I would very much appreciate it if any of you experts can confirm that I'm in the right direction w.r.t the statistical tests I need to use.   1. I'm planning to use Fisher's exact test to see if my new predictor is associated with occurrence of SCD. That is, I will divide my sample of patients into two groups, one group in which patients died of SCD, and the other in which no one died of SCD, and I will compare between the two groups the proportion of patients in which my predictor was positive. Am I using the Fisher's test appropriately here?    2. I'm planning to use McNemar's test to compare the accuracy of my new predictor with that of the conventional predictor. That is, I will have a 2x2 contingency table with the following entries:  Top left cell: # patients in which both predictors accurately determined whether or not they died of SCD. Top right cell: # patients in which my predictor was accurate but the conventional predictor was not. Bottom left cell: # patients in which my predictor was not accurate but the conventional predictor was. Bottom right cell: # patients in which both predictors were inaccurate. Is my use of McNemar's test correct?


title: Can Principal Component analyses be applied to a counting trait?

text: I am analyzing a segregating population of plants coming from an hybridization process. The experiment consists in several field plots (according to an augmented design). In each plot a segregating population coming from an hybrid plant was seeded. Therefore, the plants into each plot are segregating. I defined several traits corresponding to morphological characteristics of the plant (eg. leaves colour, flower colour, ...). The plants in each plots are segregating. Therefore the shows different characteristics for each of those morphological traits (eg. red or green leaves) and I counted the number of plants in each plot for each of those classes. Therefore I may express the data in my data-set as 'number of plants' or as percentage of e.g green/red plants on the total number of plants in each plots. Since the genetic background of the original hybrids is not known, I would like to run a PCA and a cluster analysis in order to see which populations cluster together according to those traits. Can PCA be applied to such a data-set? Which package can be used for running such an analysis in R?


title: Help in setting up and solving a transportation / traffic problem

text: Traffic light synchronization nowadays is not a tedious project, Image: ![image](http://i.stack.imgur.com/KtnTF.png) For a personal research project I'm trying to build a statistical model to solve such problems in two and higher dimensions. I would like to synchronize `x` 4 legged (approaches) intersections such that:   * The intersections are all located within the same area in a grid type format. I'd like to think of this as a 1/0 matrix, where 1 indicates a signalized intersection and 0 doesn't: **Matrix Example** [1 1 1 1] [0 1 0 1] [1 1 1 1] This assumes that all nodes are connected in a vertical or horizontal but not diagonal manner and that they could have a direction of either one or two way.   * Each Intersection has a traffic light on each approach   * The allowable movements per leg are left, right, and Through (lrt) and no U-Turns   * Each Leg can be one lane, i.e. one lane that allows for lrt or approaches can have combined, through and right (1 lane), and left (1 lane, or dedicated such as alpha left lanes, beta through lanes, and gamma right lanes. Some intersections should allow for restricted or no movement, as in an intersection could be setup as "no-left turns allowed"   * Each lane should have a distinct vehicle flow (vehicles / minutes (Time)) per movement, so a two lane (left + through, and right) approach could have the following movement: Left + Through lane: 10 vehicles / minute (through) and 5 vehicles per minute turning left Right lane: 15 vehicles per minute   * Assume green cycle time (time from beginning of green to beginning second green to depend on the vehicle volumes, yet could only be between t1 and t2 seconds. Green time = red time + yellow time). Yellow time should be assumed as a constant either 1, 2, or 3 seconds (depending on the jurisdiction). if pedestrian crossing is warranted, the red time should be assumed to be "white go ped: 4 seconds x lane" + "blinking ped: 4 seconds x lane". remember that the number of lanes would be from the opposing approaches combined.   * The system should optimize green time so that maximum flow is passed through the system at all intersections depending no the existing volumes **and** on historical values.   * Although there are historical data on the number of vehicles that enter the green network, I assume it is best to use a Monte Carlo simulation to simulate vehicle entry into the system and then calibrated to historical data.    * The end result should be an equation or system of equations that will give you the green time for each intersection depending on volumes and what to expect from the other intersections. I know this might not be the clearest question, but I'm happy to answer any questions and or comments.


title: types of baseline modelling

text: I'm searching for a model that can be used to decompose incremental effects and baseline effects and is also good at prediction (i.e. Not OLS). We have tried OLS, ARIMAX with exogenous factors, ARIMAX with exogenous factors which includes media impacts, and regional ARIMAX models (estimated regionally then summed up nationally). Maybe someone here does know some other type of model which fits our profile.


title: Comparing IRFs derived from Bayesian VARs with other extraneous information

text: I was wondering if anyone here could help me with the following: I estimate a standard Bayesian VAR with Normal-Inverse Wishart priors. I identify some policy shock in it, and then derive the IRF for variable x1 to this policy shock. This IRF for x1 is the mode of posterior distribution as suggested by Inoue and Kilian (2013). Separately, I also have information on policy surprises (e.g. like changes in the FF after FOMC meetings) and therefore I can also derive the IRF of x1 by regressing these surprises on the forward curve of x1. Thus I have a bunch of OLS coefficients that describe this same IRF for x1 though it is now identified in a completely different way. My question is how can I compare the two? Since the two models come from different probability models my idea is to perform a Wald test. For the OLS estimates I have a well defined vcov matrix for the point estimates of the IRF. No problem here. However, the IRF from my structural VAR is the modal model. To my knowledge there is no well defined notion to describe uncertainty around this mode (some sort of vcov matrix), so I have a hard understanding how the Wald test would implemented. Any suggestions?


title: Regression analysis and parameter estimates with populations

text: I've seen a little bit here about the difference between statistical inference for random samples, and what happens when we actually have population data. Most arguments seem to suggest you "never actually have the population" and the population data you think you have represents some unobservable super population which is the data-generating process. But say, we have a state, covered by counties, and we have data on some environmental aspect of each county, for example, the total area of forest. We want to see if some binary county-level outcome is related to forest - let's say presence or absence of a disease (thinking Lyme disease here). As far as I can tell, this qualifies as "population" level data. Would running a regular logistic regression, and the normal tests on parameters, yield correct estimates? One thing that comes into this is whether the "samples" are independent. In this case, forest in county A is probably related to forest in the county B next door, meaning a reduction in the degrees of freedom and imprecision in standard error estimates (as I understand it). But even if we had totally independent data, would standard, sample-based statistical inference be appropriate here?


title: Achieving statistical significance in an A/B test

text: Out of one pool, I've randomly selected a small group (group A) and put everyone else in group B. Then, I'm observing whether or not they take a desired action. If I'm looking to test at a 95% confidence, what's the best methodology and stats test to use? I read that with a large enough sample, it's ok to use a standard t-test, but can someone tell me if that's accurate? Assuming it is, I've applied it and come up with the result that group A has a much higher rate of the desired action than group B and applying the error bars, it's still showing that group A > group B (albeit a smaller rate). This is highly unlikely. At worst, the expectation was group A's actions would be more or less the same to group B. Any comments on what would cause such an unlikely result even with a high confidence level?


title: Similar sum of squares in GLM

text: I have run the General Linear Model in SPSS to analyse the effect of several demographic variables (e.g. gender, age) on the relationship between X and Y. So essentially, this is an analysis to see if there is a significant interaction between (for example) age and X on Y. (Also called moderator effect) I have found that there is a significant interaction for all the demographic variables (which makes practical sense in the context of my study). When I look at the Type III Sum of Squares output, I get the following:   * For Gender*X: 143.341   * For Age*X: 142.879   * For Education*X: 144.345 I was told that this is a bit 'unusual' because they are so close. (Unfortunately, I could not fully understand why this was unusual!) **Question 1: Why is this result unusual? What does this mean?** In addition, I have discussed that if I test the same interaction in multiple regression, the t statistics is not significant but the f statistics in General Linear Model is significant? **Question 2: Why is this the case?**


title: Logarithmic regression of form $y=a+b \log(x_1)+c\log(x_2)$ using R

text: How can I fit a logarithmic regression equation of form $y=a+b (\log (x_1)) + c(\log(x_2))$ on a data set using R? Here the main concern is that data contain zeros multiple times, so R will give infinity in the output. I have tried adding some constant to the $x$ variables, such as `log(x1+0.00001)` to avoid `Inf` but Is there any specific way to calculate this constant so that results are not affected? https://dl.dropbox.com/u/53624395/11.csv : LINK FOR DATA FILE ON WHICH I WANT TO PERFORM THE OPERATION.This is time series data and i have to perform logarithmic regression of form y=a+b(logx1)+c(logx2). and find a,b,c and then check is there any such type of relation exists or not.


title: How does one interpret the distribution over parameters in bayesian estimation?

text: I am new to Bayesian estimation. The assumption that the parameters are random variables seems a little unsettling to me. For example when considering a model for data, what physical interpretation can I provide to the equation $$ \begin{eqnarray*} P(Data) & = & \sum_{\theta} {P(Data,\theta)} \\ & = & \sum_{\theta} {P(Data|\theta)*P(\theta)} \end{eqnarray*} $$ This $P(\theta)$ i.e. probability over parameters, seems to be a bit awkward, after all how do I know what is the relative probability of the process of generation being a gaussian MM with this particular parameter combination instead of say a neural network with that parameter configuration. And further it is intuitive to think of one process generating the data, whose parameters we are guessing. But instead here we have multiple processes generating the data in tandem, i.e. a sense of a true model is lost.


title: Distribution of sum of squares of T-distributed random variables

text: I am looking at the distribution of the sum of squares of T-distributed random variables, with tail exponent $\alpha$. Where X is the r.v., the Fourier transform for $X^2$, $\mathscr{F}(t)$ gives me a solution for the square before the convolution $\mathscr{F}(t)^n$. $$\mathscr{F}(t)=\int_0^{\infty } \exp \left(i\, t\, x^2\right)\left(\frac{\left(\frac{\alpha }{\alpha +x^2}\right)^{\frac{\alpha +1}{2}} }{\sqrt{\alpha }\ B\left(\frac{\alpha }{2},\frac{1}{2}\right)}\right) \, \mathrm{d}x$$ With $\alpha=3$, the solution is possible but unwieldy and impossible to inverse to do an inverse Fourier for $\mathscr{F}(t)^n$. So the question is: has work been done on the distribution of the sample variance or standard deviation of T-distributed random variables? (It would be to the StudentT what the Chi-square is to the Gaussian). Thank you. (Possible Solution) I figured out that $X^2$ is Fisher $F(1,\alpha)$ distributed, hence will look at the sum of Fisher distributed variables. (Possible Solution) From the Characteristic Functions the average of $n-$summed $X^2$ has the same first two moments of a $F(n,\alpha)$ distribution when these exist. Hence with u the square root and doing a change of variable inside a probability distribution, the density of the standard deviation of n-sample T variables can be approximated with: $$g(u)=\frac{2 \alpha ^{\alpha /2} n^{n/2} u^{n-1} \left(\alpha +n u^2\right)^{-\frac{\alpha }{2}-\frac{n}{2}}}{B\left(\frac{n}{2},\frac{\alpha }{2}\right)}$$


title: Hidden Markov process

text: I've been reading about hidden Markov models, and I'm interested in both discrete and continuous time models (and discrete states). I have found many papers on the discrete time HMM, but not the continuous time HMM. I really want to know how to to set up the model and estimate the parameters. Does anyone know a good reference?


title: Intervention analysis with multi-dimensional time-series

text: I would like to do an intervention analysis to quantify the results of a policy decision on the sales of alcohol over time. I am fairly new to time series analysis, however, so I have some beginners questions. An examination of the literature reveals that other researchers have used ARIMA to model the time-series sales of alcohol, with a dummy variables as regressor to model the effect of the intervention. While this seems like a reasonable approach, my data set is slightly richer than those I have encoutnered in the literature. Firstly, my data set is disaggregated by beverage type (i.e. beer, wine, spirits), and then further disaggregated by geographical zone. While I could create separate ARIMA analyses for each disagregated group and then compare the results, I suspect there is a better approach here. Could anyone more familiar with multi-dimensional time-series data provide some poitners or suggestions?


title: Popularity of Bayesian methods in statistics and machine learning

text: I know I have seen some research, perhaps in the contexts of time-varying topic models, on the popularity of Bayesian methods in statistics and machine learning over the last 20 years. Unfortunately I cannot seem to find the citation. Does anybody know a paper that includes a plot over time of the popularity of Bayesian methods in published statistics and machine learning journals?


title: Reference Request: Generalized Linear Models

text: I am looking for an introductory to intermediate level book on Generalized Linear Models. Ideally, in addition to the theory behind the models, I would want it to include applications and examples in R or another programming language - I hear SAS is also a popular choice. I intend to study it on my own and so it would help if it provided the answers to its own exercises. You can assume I have taken the traditional year-long courses in calculus and probability theory. I am also familiar with the basics of regression analysis.


title: Combining probabilities for sports prediction

text: I use two criterion ($w$ and $x$) for predicting the outcome of football matches. Analysis of historical records has provided me with two best fit linear equations. The probability that the team playing at their home ground will win based on criteria w, is $y_1 = 0.02w + 0.38$. A different probability based on $x$ is $y_2 = 0.01x + 0.37$. What method should be used to combine these two equations to produce a single equation to obtain $y$, or should an entirely different approach be used to arrive at a single equation?


title: Data mining algorithm suggestion

text: I would like to use data mining to try to find a good workout schemes. The input dataset will contain the parameters of a set of workouts with dates and different performance and medical measures. The problem is that the influence of each individual workout will be different depending on the time that has passed after this workout. For instance on the next day after the workout the performance will degrade, but eventually will improve to the higher level than before the workout. So each performance measure will be the cumulative result of different workouts. I know that I can run a series of workouts of the same type during a period of time and then use difference in the performance as a predictor. But I was wondering if there are any algorithm that allow to take into account the time component and run analysis that will take into account individual workouts.


title: Choosing a distribution to deal with over-dispersion

text: I'm working on a life insurance problem: trying to simulate the total dollar amount of claims in a year. To do this, I have a record for each person that contains their amount of insurance and an estimate of the probability that they will make a claim during the next year. I have run thousands of one-year simulations and developed the range of outcomes. Unfortunately, the real historical experience is far more volatile than my analysis suggests. Almost every year looks like a 1-in-100 year. I used the binomial distribution (thinking simplistically that either a person does or doesn't make a claim), but I have been advised that my data suffers from over- dispersion. That sounds possible, because obviously we do not really know each person's precise probability of making a claim. Our estimate is wrong for any particular individual, but reasonably accurate for each sub-group. The recommended solution was to use a negative binomial distribution so that the variance can be calibrated separately from the mean. Certainly this property would be helpful, but literature on the negative binomial distribution focuses on applications of counting successes and failures. That doesn't seem relevant to what I'm trying to do. Is this really a good probability distribution to use in this context, and why?


title: Fit GAMs with minimum constraints

text: I have a two-dimensional predictor plane, 0-1-Observations and a priori knowledge of minimal probabilities per combination of the predictors. I would like to fit a model (e.g. GAM from the mgcv package) that yields predictions that lie above the minimal values. The problem is how to specify the model. The following is a toy example where the minimal probability is 0.5.               rm(list=ls())     require(mgcv)     set.seed(1)     l<-1000     x<-seq(0,1,l=l)     p<-0.5     y<-as.numeric(runif(length(x))<1-x^4)          df<-data.frame(x=x,y=y,min=p)     m1<-gam(y~s(x)+min,family=binomial(link="logit"),data=df)     predict(m1,newdata=data.frame(x=1,min=0.5),type="response")      I found that offset does not help. Subtracting the minimal probability from the response does not help as this yields negative responses which do not make sense. How can I adapt the model to accomodate for the given minimal value?


title: What is simply meant by reduced form?

text: In econometrics, what is meant by reduced form? Also, what are people looking for when they say "I would like to see the reduced form estimates." This has been thrown around at work and individual explanations and Google searches are overly technical. Hoping someone where would be able to give a simple example.


title: Compare errors of two methods

text: I wish to compare two regression methods, which has been run 20 times and each time both methods are given the same input (paired). For example I might observe the errors: > Method A: 1.2 0.9 1.3 1.5 1.2 > > Method B: 1.3 1.1 1.4 1.5 1.1 How do I test which method is better? When I say better first thing that comes to mind is test $\mu_A<\mu_B$. Are there other ways of testing what is 'better'. For the moment I am doing a paired T-test.


title: How to ensure data is normally distributed for the purpose of performing a continuous wavelet transform?

text: I have a vector which contains several values e.g.               a <- runif(8760)      I wish to perform continuous wavelet transform to the data but the data is required or preferred to be normally distributed. Can anyone provide any information regarding how I could transform the data to be normally distributed. With my own data set I have tried to standardise the data and perform a log transformation but these dont seem to help. I have tried:               b <- log(a)     c <- scale(a)


title: What is the error term in a Multiplicative Error Model?

text: In the Multiplicative Error Model (MEM) specification by Engle in "New Frontiers for ARCH models", he wrote that MEM specified an error that is multipled times the mean and the specification is: $$ x_t = \mu_t \epsilon_t $$ $$ \epsilon_t | \Im_{t-1} \sim D(1,\phi_t^2) $$ I am not sure I understand what epsilon is supposed to be, is it a normal of mean 1, but then what is the variance supposed to be ? Engle mentions that the MEM is good to model the process which have to be positive. From my understanding this is only true if $\epsilon_t$ is positive, so unless $D$ is a log-normal distribution or another distribution guaranteeing positivity, you cannot guarantee that $x_t$ itself will be positive, is my understanding correct ?


title: Time series classification pointers

text: I'm new to ML/statistical learning and would like a few pointers in what I need to study to solve my problem. (FTR, I've done an Intro to AI course, classifying Fisher's Iris dataset and things like that). I'm trying to classify regions or slices in time series data. These are associated with a label. I have many examples of time series with slices and labels (i.e a training set). For new data, I would like to associate time slices with a label. However, a label for a certain time slice will also depend on previous estimated labels in the same set. Forgive the vague terminology, but I would like to be pointed in a particular direction for a problem of this type. I've heard Hidden Markov Models and Dynamical Systems would be relevant to study for this task.


title: How to interpret  two standard deviations below the mean of a count variable being less than zero?

text: I asked people how many times they visited their local pub in a 'normal' week. The result can be zero, one, two, three, four, and five and more.   * The mean is 2 and the standard deviation is 1.3.    * So two standard deviations above the mean is 4.6.    * However, two standard deviation below the mean is -0.6. Is this negative figure an error? How do I interpret it?


title: Logistic regression with partially observed variables

text: Let's say I have a binary variable to explain using "some" logistic regression approach. The set of variables I have at disposal are exhaustive for a given period of time and for a longer period a subset of those data are missing. I think this is a very common issue and was looking for some references about the classic (or maybe more exotic) methodologies available to deal with this kind of situation. Ideally, I would like to fit the model to all the available data without filling the missing value with best estimates. I was thinking that maybe nesting model or conditionning appriopriately would be possible.


title: unbalanced haar wavelets

text: My aim is denoising of the smooth functions based on second generationwavelets.In (unbalanced haar technique for nonparametric function estimation) signal are denoised using unbalanced haar wavelets. I coud understand best top-down unbalanced haar decomposition,but it no use for the aim of smooth functions.and it suggested by the author of the article ‘’bottom-up’’technique for smooth function,could some one please help me with ‘’bottom-up’’ technique? And also could some one please help me with ‘’design- adapted wavelet’’described in ‘’smooth desined adapted wavelets for nonparametric stochastic regression’’.in this paper with the tool of lifting predict and update filters has been designed .so has changed the vanishing moments.this prt has utilized an algorithm than I have no idea about. Help are appre ciated


title: Avoiding large variances when taking the logs of small values

text: I have two random variables $(X$ and $Y)$ that are always positive. The assumption I'm making is that their logs follow normal distributions (i.e., $N(\overline{\log(X)},s^2_{\log(X)})$ and $N(\overline{\log(Y)},s^2_{\log(Y)})$) . I'm interested in the difference between them, and therefore I'm assuming that the difference of their logs is also normally distributed $(\log(X)-\log(Y) = \log(X/Y) \sim N(\overline{\log(X/Y)},s^2_{\log(X/Y)})$). My problem is that $X$ and $Y$ can take values that are very close to $0$ (as small as $10^{-30})$, and in that case the variances of their logs become extremely large. For example, if the mode of $X$ is around $10^{-4}$, but $X$ has a left tail down to $10^{-30}$ and a right tail up to 1.5, its log will have a huge variance and obviously also that of $\log(X)-\log(Y)$. My questions is how to deal with this. Probably something in the lines of adding a constant to $X$ and $Y$, but I'm not sure what that constant should be. I thought of adding an $\epsilon$, but that won't help the case of the example of $X$ I mention above.


title: Linear Regression Intuition behind least squares

text: I try to fill the gap of linear regression by following the textbook **Patterns Recognition and Machine Learning by Cristopher M. Bishop**. It's known that the error function for linear regression is $E(w)=\frac{1}{2} \sum_{n=1}^{N} {y(x_n,w) - t_n}^2$ (defined on page 5) where $y(x_n,w)$ is a hypothesis, a linear function with respect to $w$. My question is why this a right function to minimize, and the intuition is given on page 28, however I found it very vague. > ...we shall assume that, given the value of $x$, the corresponding value of > $t$ has a Gaussian distribution with mean equal to the value $y(x,w)$ of the > polynomial curve given by (1.1). Thus we have $p(t|x,w,\beta) = > N(t|y(x,w),\beta^{-1})$ where $\beta$ is a inverse variance of the > distribution. The further explanation is straightforward, The question is why exactly the hypothesis $y(x,w)$ we take as a mean, for me it's not intuitive why we can make such an assumption. I will be very helpful if could explain it on the basic level. In my opinion the key point to the understanding is the following figure, which is attached to the explanation, unfortunately I still don't understand what it depicts. ![enter image description here](http://i.stack.imgur.com/mPXle.png)


title: Are there models for "censored" spatial point processes?

text: This is a problem I'm encountering in the context of analyzing a data set comprised of all crime locations in a city over a fixed time interval, although it could potentially arise in other types of point processes. The problem has to do with the fact that crime locations are not observed _exactly_. They are located by whatever street address the police wrote down on their report. In my city of study (and probably almost any city), there are "natural barriers" in the landscape that make these "locations" inherently imprecise. For example, suppose there is a large nature reserve/park in the middle of a city. Then, all crimes occurring within the nature reserve are mapped to a single address - presumably the address of the main office of the nature reserve. This kind of "censoring" causes artificial clustering in the data set and, most likely, biases the estimate of the intensity function and the associated covariate effects, etc. I'd imagine these sort of "natural barriers" exist in pretty much any city so that this issue has probably been fomented by other researchers so my question is: **Are there known methods for handling this type of data?** At the moment, I have analyzed this data set using an ordinary inhomogeneous poisson process model and have gotten some interesting results. I really think these results are "real" based on previous descriptive analysis, etc. and the model fits reasonably well except near these natural barriers where there are **huge** residuals due to smoothing over "zero density" areas, making the model wildly fail any kind of "goodness of fit" test despite the fact that the empirical density of the observed data agrees fairly closely with the empirical density of data simulated under the fitted model. Here are the main possibilities I've considered (and why I've decided against them):   * Delete these "natural barriers" from the window of observation and view this as a point process on a grid. I don't want to do this since it fundamentally changes the parameter you're estimating and, effectively, sweeps the "censoring" issue under the rug.    * Bin the data into areal units (e.g. census-based groups), since this censoring infrequently crosses over census boundaries. This may be a good solution in some cases, but oversmoothing is a concern and, more importantly, the city I'm working with is too small (not enough census units) to do this.    * Develop my own model for this. I'm pretty sure this is what I have to do but I wanted to check first that I'm not reinventing the wheel. Based on my own literature search, I'm not, but some specialists here may know something I don't..


title: Visualizing large dataset with multiple subgroups

text: I have a large data frame in the following form (I apologize for this formatting):               Site    Season  T          SC    pH    Chl   DO.S   DO      BGA  Tur    fDOM    Flow    Rainfall    Solar      Rain     300N    Winter  14.05   1692.77 7.93    NA  82.26   8.42    NA  9.25    NA      NA      0.00          219.18     no      If you can't understand the formatting, there are 12 numerical factors, and 3 categorical factors (`Site`, `Season`, `Rain` [yes/no]). Each row represents the average daily values that I have calculated from 15-minute time series. I have spent a good amount of time doing data exploration (linear regression analysis, looking at time series plots for patterns), but haven't found a method that works for me yet. I have also worked with `corrplot`, correlation matrices, and covariance functions in an arduous way, where I subset each categorical combination and found `corrplot`s for each (I have also tried it with `ddply`, but the resulting format is not in the correlation matrix format that is easy to plot). I have also attempted PCA on the data to little avail. My question is first and foremost, does anyone have an idea for data visualization of this kind of dataset? The main question I am after is, "What are the factors that influence `DO` (dissolved oxygen)?". How does this change by location (`Site`), `Season`, and with the influence of `Rain`. I would really like a quick method for shooting out correlation matrices (or heat maps; I have tried both) for each categorical subset. I tried this with `ggplot` and `facet_wrap`, but it wasn't happening for me. I also tried `ggpairs` from the GGally package, but honestly didn't spend too much time with that method. I was starting to get into the idea of star graphs (on polar coordinates), which can be used to visualize repeating periodicity in time series, but am running out of time and decided to seek the advisement of Stack Overflow. I really appreciate any advice or thoughts on visualizing this data that come to your mind. I feel like some combination of `ddply` and graphing is what I need, but I haven't gotten there yet. Thank you for your time. EDIT: `dput` of the data frame in question:               structure(list(Site = structure(c(2L, 2L, 2L, 2L, 2L, 2L), .Label = c("2100S",      "300N", "3300S", "800S", "Burnham", "Center"), class = "factor"),          Season = structure(c(4L, 4L, 4L, 4L, 2L, 2L), .Label = c("Fall",          "Spring", "Summer", "Winter"), class = "factor"), T = c(14.05,          14.18, 14.5, 14.58, 14.07, 11.91), SC = c(1692.77, 1671.31,          1680.71, 1661.79, 1549.56, 1039.63), pH = c(7.93, 7.92, 7.96,          7.95, 7.93, 7.79), Chl = c(NA_real_, NA_real_, NA_real_,          NA_real_, NA_real_, NA_real_), DO.S = c(82.26, 78.79, 82.05,          80.92, 74.33, 73.96), DO = c(8.42, 8.04, 8.31, 8.18, 7.61,          7.97), BGA = c(NA_real_, NA_real_, NA_real_, NA_real_, NA_real_,          NA_real_), Tur = c(9.25, 9.77, 9.41, 10.6, 40.38, 50.25),          fDOM = c(NA_real_, NA_real_, NA_real_, NA_real_, NA_real_,          NA_real_), Flow = c(NA, 178.08, 178.53, 188.13, 306.15, 382.22         ), Rainfall = c(0, 0, 0, 0, 0.01, 0.81), Solar = c(219.18,          228.33, 244.3, 247.69, 105.15, 220.73), Rain = structure(c(1L,          1L, 1L, 1L, 2L, 2L), .Label = c("no", "yes"), class = "factor")), .Names = c("Site",      "Season", "T", "SC", "pH", "Chl", "DO.S", "DO", "BGA", "Tur",      "fDOM", "Flow", "Rainfall", "Solar", "Rain"), row.names = c(NA,      6L), class = "data.frame")


title: How to count expected number of paterns?

text: Member of parliament can vote or can be absent during voting. Let assume: If, for a given MP, there is a pattern present-absent-present (`010`) in consecutive voting, this absence is strategic - MP does not want to reveal his preferences. For each MPs I have a data in the format:                   00011001000101_00001000_10_011_00000      Where `0` means presence of MP during voting, `1` - absence. `_` means there was a interruption (ie. coffee break, night etc.) between voting. I can count probability of absence of MP - it is number of all absences divided by a number of all votings. Assuming absences are independent, and probability for one absence is known (computed as above), how can I count expected number of `010` patterns in data? For above data, I want to do it this way:                   p=9/(23+9)=0.28125 - probability of absence         (1-p)*p*(1-p) = 0.1453 - probability of 010 pattern         number of all possible length 3 patterns in given data:         (14-2)+(8-2)+(2-2)+(3-2)+(5-2) = 22         so the expected number of 010 pattern is 22*0.1453 == 3.2      is my reasoning OK? How can I test - based on expected and real numbers of `010` patterns - if this is only by chance or not?


title: How to calculate max/min scales on a scatter plot

text: I have 3 log scatter plots that I want to establish smooth maximum and minimum lines. What is the usual mathematical method to do that? (Image and Excel file links below.) The black lines on the scatter plot images are hand drawn. The third scatter plot is especially tricky and not amenable to a moving average plus stddev because of the data clumping. Note: This is time series data so new data constantly comes in. In other words, I cannot just use the whole data population in one shot. Any ideas would be greatly appreciated. Excel File: https://dl.dropboxusercontent.com/u/44057708/Three%20Scatters.xls Image at: https://dl.dropboxusercontent.com/u/44057708/ThreeScatters.jpg ![enter image description here](http://i.stack.imgur.com/lxHc2.jpg)


title: Sub sampling a dataset without losing the distribution

text: Suppose that I have a huge data set with high number of dimensions (components). For computational purposes I would like to take just a sample of this data set and work with it instead of working with the huge data set itself. Is there a good way or technique to sample a dataset? I was thinking to choose a sample randomly but in this way it is possible that the distribution of the sample will not be the same as the distribution of the original data set. Example: If I have 2000000 data on that data set and within these data 80% of those data will be from one group(cluster) and the others will be in the other clusters (not related to the 80% of above mentioned data). Now, if I randomly select a sample most probably I will have all almost the data taken from the 80% and not from the other 20% of the data. So I am looking for a good way to take a sample from such a data set and I dont want the small clusters to be lost. Because as I said those 80% of original data belong to one closter and the remaining 20% belong to some other small clusters. As a result, I think that if I randomly take a sample from these data and then I apply a clustering algorithm on this sample, most probably I will be loosing small clusters