title: Average value from ranges in a survey

text: I have a questionnaire where some of the questions concern ranges, the answer alternatives are:   * 0   * 1 - 15   * 16 - 30   * 31 - 45   * 46 - 60   * 60+ (the numbers are minutes per week, not that it matters) If I would like to present the result as an average number, how would I do this? My thinking right now is just to take the average in the range (1 - 15 => 7), times how often it was selected. Then sum it for all alternatives and divide by the total number of answers. Not sure how to handle the 60+ value yet though. But it seems to simplistic and since it looks like a problem that must occur quite often, there might be some better method. **Edit:** A little background: The survey is regarding sources of inefficiencies in the workplace, like issues with slow servers, slow computers, using the wrong programs, etc. There are questions for several areas, and for each area the user will answer an estimate of how much time is lost, in minutes per week, in that particular area. I was looking for an average since I thought it would be an effective way to present the results; "On average, a developer loses X minutes per week due to an inefficient workplace. This means we could save Y money if we....". There is no need for it to be precise, it will just be used as an indicator and give some incentive for change.


title: Are mean normalization and feature scaling needed for k-means clustering?

text: What are the best (recommended) pre-processing steps before performing k-means?


title: Joint distribution of two gamma random variables

text: I am so puzzled by this problem. Given two variables $X_1$ and $X_2$, such that $X_i \sim \mathrm{Gam}(a_i, b)$, find the joint distribution of $X_1$ and $X_2$. I understand how to proceed if variables are independent, but the general case is a bit unclear. Thanks for your help!


title: Comparing non-nested models with out of sample likelihood

text: I recently read a paper in which the authors claim that in order to compare the forecasting performance of two non-nested models, models A and B, a valid procedure is to fit models A and B on the same data set, and compare the average likelihood of the fitted models computed on the hold-out data set. All that is required is that models A and B are expressed as probability densities for the same variable. I am using the language of the paper: strictly speaking the quantities being compared are fitted models evaluated at data points, not likelihoods. Accepting this likelihood as the forecasting accuracy metric, this method for checking forecast accuracy has some intuitive appeal, although I have misgivings. Can these likelihoods be meaningfully compared without some sort of normalization? The out of sample probabilities will not add up to one. I'm not able to come up with a straightforward example in which this test would give a spurious result. update: I was able to produce a simple example in which direct comparison of out of sample likelihoods gave a misleading result: let the dependent variable y be a linear trend plus normally distributed error term, and let the explanatory variable x be a linear trend plus independent normally distributed error term. Generate 100 points for y and x. Model A is a linear regression, while for model B I used a regression with Student t-distributed errors (with two degrees of freedom). I trained on the first 50 points and tested on the second set of 50 points, and repeated with training and test sets interchanged. I repeated with three choices of variance for the data generating process. Model B gave higher average out of sample likelihoods in all cases. This example is a bit contrived but does illustrate my concern.


title: Which test to use to compare means/medians when one variable has normal distribution and the other does not?

text: Which 2-tailed test is best to use to compare means/medians when one variable has a normal distribution and the other does not?


title: Difference between ANCOVA and Hierarchical Regression

text: Is there a difference between ANCOVA (as performed under the 'General Linear Model (GLM)') and Hierarchical Regression (as performed under 'Regression') in SPSS? I am testing the main effects and interaction of X1 (continuous) and X2 (categorical) on Y (continuous). I understand GLM incorporates regression and one of the main advantages of using GLM over the regression function is that GLM (factor) creates dummy coding. In regression, this must be done prior to inputting the categorical factor in the analysis box. The same also applies to the interaction product. I personally find GLM (ANCOVA) output easier to understand.


title: Cross validation and parameter tuning

text: Can anyone tell me what exactly a cross-validation analysis gives as result? Is it just the average accuracy or does it give any model with parameters tuned? Because, I heard somewhere that cross-validation is used for parameter tuning.


title: Regression specification choices

text: I am studying a population of individuals who all begin with a measureable score of interest (ranging from -2 to 2) [call it "old"], then they all undergo a change to a new score (also ranging from -2 to 2) ["new"]. Thus all the variation is in the change (which can be positive or negative), and there are also a variety of predictors that help to explain variation in the amount of change. My initial model is simply:               change = a + bx + e      where x is my vector of predictors. But now I'm concerned that some of these predictors could be correlated with the baseline (old) score. Is this, then, a better specification?               change = a + bx + old +  e      Or perhaps               new = a + bx + old + e      Thanks!


title: Might be an unbalanced within subjects repeated measures?

text: I ran a within subjects repeated measures experiment, where the independent variable had 3 levels. The dependent variable is a measure of correctness and is recorded as either correct / incorrect. Time taken to provide an answer was also recorded. A within subjects repeated measures ANOVA is used to establish whether there is significant differences in correctness (DV) between the 3 levels of the IV, there is significant. Now, I'd like to analyze whether there is significant differences in the time taken to provide the answers when the answers are 1) correct, and 2) incorrect. My problem is: Across the levels there are different numbers of correct / incorrect answers, e.g. level 1 has 67 correct answers, level 2 has 30, level 3 has 25. How can I compare the time take taken for all correct answers across the 3 levels? I think this means its unbalanced? Can I do 3 one way ANOVAS to do a pairwise comparison, while adjusting p downwards to account for each comparison? Thanks


title: How to utilize the result of aov function?

text: Based on the responses from the question How to draw a side-by-side plot mentioned in "Graphical Display as an Aid to Analysis", I am considering implementing the function myself. The problem is that I don't know how to utilize the result of `aov()` function. I have to see how many factors there are, the names of factors, what are their main effects and interaction effects, etc. I can get residuals from `aov(y~x)$residuals`. I consulted the documentation for `aov()`, but I still have no clue how to do this. There must be some other way much complex than something like `t.test(x,y)$p.value`. Would you help me?


title: What is the expected number of runs of same color in a standard deck of cards?

text: Standard deck has 52 cards, 26 Red and 26 Black. A run is a maximum contiguous block of cards, which has the same color. Eg.   * (R,B,R,B,...,R,B) has 52 runs.   * (R,R,R,...,R,B,B,B,...,B) has 2 runs. What is the expected number of runs in a shuffled deck of cards?


title: How to simulate a custom power analysis of an lm model (using R)

text: Following to the recent questions we had here. I was hopping to know if anyone had come across or can share **R code for performing a custom power analysis based on simulation for a linear model?** Later I would obviously like to extend it to more complex models, but lm seems to right place to start. Thanks.


title: Monomial distribution of $X^a \cdot Y^b$

text: What is the distribution of the following monomial? $$X^a \cdot Y^b$$ where $X$ and $Y$ are normal random variables and $a$ and $b$ are natural numbers. For example, when $X \sim N(0,1)$, $a=2$, and $b=0$ it is a Chi-squared distribution, which has a variance of 2. What if we have $n$ independent variables $X_1, X_2, \dots , X_n$, with $X_i \sim N(0,\sigma^2)$ and some natural numbers $p_1, p_2, \dots,p_n$. What can we say about the variance of the following r.v.? $$X_1^{p_1} \cdot X_2^{p_2} \cdots X_n^{p_n}$$


title: Statistical power of unit root test

text: When studying time series, I once heard the statement that > unit root test is less powerful. I hereby have two questions:   1. What does it mean for a test to be powerful?   2. What causes the unit root test to be less powerful?


title: How useful is the CLT in applications?

text: My lecturer just covered the **Lindberg-Levy central limit theorem** and the multivariate version, the **Lindberg-fuller CLT**. I understood the basic concept and I can derive it, etc. But it would help my understanding a lot if someone could explain how all this used in real life applications of **econometric** analysis? Ive read some claims that the CLT is nice only on a piece of paper. Some really cool industry applications or references would be appreciated.


title: What is the relation between Kullback loss and L1 and L2 loss?

text: I tried to find some relation between these distance (loss) measures, but couldn't find any references. However, I think it must something like this: $$ \sqrt 2*D_{KL} < L_1 < L_2 $$ Is that right?


title: Comparing time series: Pearson correlation, Kendall's tau b or Spearman's rho?

text: Which correlation coefficient is the most appropriate to compare 2 time series? I want to compare the variation of one variable for 2 regions, have regional data for the last 30 years. Is Pearson correlation ok or should I rely on Kendall's tau b or Spearman's rho and why? I tried to google it and analyse what I found, but I'm still not sure.


title: AIC for semi-parametric model

text: I have learnt the definition of AIC for parametric model. But, what I want to know is that is there a semi-parametric version of AIC? Have you heard about it? And if you know, please show me a link, so I can have a look.


title: Is there a simple method to calculate the minimum number of datalines needed for a machine learning algorithm with x variables?

text: Here is the problem: I am making a machine learning algorithm that takes the inputs and outputs of some software I've written, and I don't know how many datalines to produce to get results that are a 'good' fit. I realize the answer is 'the more the better' but I'm looking for any sort of minimum requirement. I also realize that the greater the number of variables the greater the number of datalines required. So I'm looking for a rule of thumb for number of variables to number of datalines.


title: How to implement a weighted 3-way ANOVA in unbalanced design?

text: I have implemented a three-way anova with type III sum of squares in c++. Since some of my experiments (observations) are more important (more informative), I want to give them a higher weight in my analysis. For example, an experiment which is very important has a weight of 10, and a relatively important one has a weight of 5, and so on. To implement it, I repeat such observations in respect to its corresponding weight, 10 times, 5 times, ...   I used the same concept in 2way anova, but there I used the conventional formulas of Sum of Squares, because my design was balanced. So there I just multiplied the value of each observation and the number of times I have seen that to the weight.   Here, repeating the items makes the design matrix so big and increase the computational complexity. Now the problem is that what if I don't want to repeat them, but use a weight matrix? and what if the weights are not integer values (so I really can't repeat one item 0.3 times)?   I found this formula here:               H^=X(X'WX)-1XW      So I put my weights into W matrix and used this formula. To check if it works, I used the conventional method and repeated observations as often as the weight value, and gave it to MATLAB. But I got different SS values. Could you please kindly tell me how I should change this formula? Thanks,


title: Repeated 100x 10-fold cross validation, what is the sample size when doing an significance test?

text: I iterated my 10-fold cross validation 100 times for several methods. Now I want to use a t-test to test if the results are significant. However, I'm not sure what the sample size is. Is the sample size the original amount of samples, or is it the original amount of samples x 100? ## edit: For university we need to classify 3 cancer types and give an estimation of how well our model will perform. We received a dataset with 100 samples. We split the data up into a training and test set using stratified sampling with a ratio of 0.3 and 0.7. The resulting training set consists of 69 samples, and the test set out of 31 samples. We used repeated cross validation because of this paper: http://www.cse.iitb.ac.in/~tarung/smt/papers_ppt/ency-cross-validation.pdf The repeated cross-validation is done on the same training set, but with the folds are randomly chosen every time, so they should be different every time. The significance we want to test is if the accuracy of one model is significantly better than the accuracy of a different model.


title: (Prediction) constrained OLS in R

text: I have a collection of $n$ datapoints $(y_i,\bf{x}_i)$ in $\mathbb{R}^{p+1}$ and would like to estimate the following model in `R`: $$\underset{\bf{b}\in\mathbb{R}^{p}}{\arg.\min}\;\sum_{i=1}^n(y_i-\bf{x}_i'\bf{b})^2$$ $$u.c.\;\;\;0\leq \bf{x}_i'\bf{b}\leq 1\;\;\;\forall i$$ anybody has a pointer to an efficient way of doing this? is there a re- parametrization of the OLS problem that would allow for this?


title: Comparing effects y1 ON x vs. y2 ON x in different regressions

text: Let's assume we have predictors x1-x5 and dependent variables y1-y9 in one dataset. We have a certain hypothesis about x1: it should have differently strong effects on our 9 dependent variables y1-y9. We perform 9 regressions, and find that x1 significantly predicts all y with p<.001. Now we want to find out whether these (highly significant) effects are different from each other (just because they are significant effects does not mean they are equally strong). Two questions: (1) What information that a regression provides us with would give us insight into this question? Unstandardized beta? Standardized beta? t? (2) Is there a statistical test we can perform to find out whether the strengths of predictions are different between (y1 ON x1) and (y2 ON x1) and (y3 ON x1)? Software I can use is R and SPSS; field is psychology/medicine.


title: Program Impulse Response Functions for VAR

text: I'm trying to program impulse response functions for a VAR model using Cholesky decomposition. The thing is I do not completely understand how I should do this when I read in the literature. Suppose I have: $$ \begin{bmatrix}x_t\\y_t\\z_t\end{bmatrix}=\begin{bmatrix}\alpha_1\\\alpha_2\\\alpha_3\end{bmatrix}+\begin{bmatrix}b_{11}&b_{12}&b_{13}\\b_{21}&b_{22}&b_{23}\\b_{31}&b_{32}&b_{33}\end{bmatrix}\begin{bmatrix}x_{t-1}\\y_{t-1}\\z_{t-1}\end{bmatrix}+\begin{bmatrix}u_{1t}\\ u_{2t}\\ u_{3t}\end{bmatrix} $$ which we can write as $$ \mathbf{x}_t=\mathbf{a}+\mathbf{B}\mathbf{x}_{t-1}+\mathbf{u}_t. $$ Further, suppose the covariance matrix of $\mathbf{u}_t$ is $$ Cov(\mathbf{u}_t)=\Sigma_u=PP^\prime. $$ Now, let's say I want the impulse responses to a unit shock in $u_{1t}$. I want the effects then on say $\mathbf{x}_t, \, \mathbf{x}_{t+1}, \, \mathbf{x}_{t+2}$ and $\mathbf{x}_{t+3}$. As I understand it, the orthogonalization is done by multiplying the error vector with $P$. Let's call the responses in period $t$ to the shock $\mathbf{x}_t^*$. Would what I am interested in then be (assume unit shock in $u_{1t}$ such that $\mathbf{u}_t^*=\begin{bmatrix}1&0&0\end{bmatrix}^\prime$): $$ \mathbf{x}_t^*=P\mathbf{u}^*_t=P\begin{bmatrix}1\\0\\0\end{bmatrix}\\ \mathbf{x}_{t+1}^*=\mathbf{B}\mathbf{x}^*_t=\mathbf{B}P\mathbf{u}^*_t\\ \mathbf{x}_{t+2}^*=\mathbf{B}\mathbf{x}^*_{t+1}=\mathbf{B}\mathbf{B}P\mathbf{u}^*_t\\ \mathbf{x}_{t+3}^*=\mathbf{B}\mathbf{x}^*_{t+2}=\mathbf{B}\mathbf{B}\mathbf{B}P\mathbf{u}^*_t $$ Now, extend this to include more lags (for example 4). The model is then $$ \mathbf{x}_t=\mathbf{a}+\sum_{k=1}^4\mathbf{B}_k\mathbf{x}_{t-k}+\mathbf{u}_t. $$ Thus the impulses are: $$ \mathbf{x}_t^*=P\mathbf{u}^*_t=P\begin{bmatrix}1\\0\\0\end{bmatrix}\\ \mathbf{x}_{t+1}^*=\mathbf{B}_1\mathbf{x}^*_t=\mathbf{B}_1P\mathbf{u}^*_t\\ \mathbf{x}_{t+2}^*=\mathbf{B}_1\mathbf{x}^*_{t+1}+\mathbf{B}_2\mathbf{x}^*_t=\mathbf{B}_1\mathbf{B}_1P\mathbf{u}^*_t + \mathbf{B}_2P\mathbf{u}^*_t\\ \mathbf{x}_{t+3}^*=\mathbf{B}_1\mathbf{x}^*_{t+2}+\mathbf{B}_2\mathbf{x}^*_{t+1}+\mathbf{B}_3\mathbf{x}^*_{t}=\mathbf{B}_1\mathbf{B}_1\mathbf{B}_1P\mathbf{u}^*_t + \mathbf{B}_1\mathbf{B}_2P\mathbf{u}^*_t+\mathbf{B}_2\mathbf{B}_1P\mathbf{u}^*_t+\mathbf{B}_3P\mathbf{u}^*_t $$ Is this line of thinking correct? If so, then this simple R code should be fine:               library(vars)     set.seed(1)     x <- rnorm(100)     set.seed(2)     y <- rnorm(100)     set.seed(3)     z <- rnorm(100)     data <- cbind(x, y, z)          model <- VAR(data, p=4, type = "const")     u <- matrix(c(1, 0, 0), ncol=1)     P <- chol(cov(residuals(model)))     B1 <- Acoef(model)[[1]]     B2 <- Acoef(model)[[2]]     B3 <- Acoef(model)[[3]]     B4 <- Acoef(model)[[4]]          xt <- P %*% u     xt1 <- B1 %*% xt     xt2 <- B1 %*% xt1 + B2 %*% xt     xt3 <- B1 %*% xt2 + B2 %*% xt1 + B1 %*% xt      Any input would be very much appreciated!


title: How do you predict the value of new instance, when the training data were normalized?

text: I estimated a Partial Least Squares model where the X matrix had normalized columns. Now I want to predict the value for a new instance (which is a frequency vector summing to one.) I assume that if I just use the raw frequency values, the predicted value won't be on the same scale as the scenario where my 'new' instance was taken from the normalized X matrix. (i.e. Comparing the fitted values of the model with predicted value of new instance.) I was thinking of adding the new instance as the bottom row of the original non-normalized X matrix, normalizing, and then using the values from this new bottom row to predict. Alternatively, I could standardize by using the column means and standard deviations from the original non-normalized X. Is one method preferred to the other? Is there a better way?


title: Learning an extra programming language (other than R)?

text: My work involves building statistical / econometric models using R, SPSS modeler. I am also doing my PhD (part time) in econometrics. In order to do more advanced data / model visualisation I am thinking about to pick up another programming language. Any suggestion will be much appreciated.


title: Pre and Post, treated and un treated but from same subject

text: A treatment was given to one hand of a subject, and a single outcome metric is measured for both hands, twice pre and several times post treatment. What is best practice for assessing effectiveness of treatment? Treated and Untreated "groups" really are paired.


title: Designed experiment where tool wear is an important factor

text: I have to design an experiment where tool wear is a factor (not a response). I'm not trying to study tool wear minimizing. I'm trying to study the effect of tool wear along with 3 other factors on the number of a certain defect type. What kind of design should I use? I considered dividing the tool wear into 3 conditions (low, med, high) and then running each treatment within the conditions. But I'm not sure that is a correct method. Tool wear is considered to be a highly significant contributor to the defect rate. Interactions are also important here. This is a forging operation. The other factors are billet temp, spotting location and pancake height. I have plenty of replication opportunities and plan at least 10 parts for each treatment. Many Thanks - Kelly


title: A "Gambler's Loss function"?

text: What is a good loss function for a predictive model used by gamblers? I've been reading a bit about loss functions recently. I've always just went with MSE (e.g., for a couple of neural network projects) and didn't ask questions. I didn't realize exactly how arbitrary MSE actually is. And speaking of this, could somebody explain a simple practical situation where MSE can be derived as the "correct" loss function? Anyway, I came across Log-Loss (which is the same as cross entropy) and am interested because I'm curious about some probabilistic models. I understand how one derives this loss function from information theory, but when we talk about probabilities, we're often involved in some form of gambling on an outcome. It's not clear to me that efficient transmission of information translates to the type of utility one is usually looking for in a predictive/probabilistic model. If, for example, I had a model that was meant to predict the winner of the 2012 US presidential election, and I had used this to move money around on a future's market like Intrade, how might I determine a loss function for my prediction - assuming I'm able to continue making bets as the market fluctuates? Same type of thing should apply for any market where I am able to make handicapped bets on the occurrence of some event. Or is this really regret, and is regret completely different than a loss function?


title: Comparing pre-retrofit utility bill time series data to post retrofit data

text: I have hourly utility bills for two years. The first year's data are from before a retrofit was performed on a home and the second year is after the retrofit was performed. How do I compare the two time series statistically to make claims of demand savings during a specific time period of the day. Can this be done in R, if so how?


title: Flipping random coins from a bag - equivalent to a single coin?

text: My first and I think naive question here. I am trying to model a certain business, and the simplest model I am willing to test is: 1. there is a bag of differently biased coins. 2. every step, a single coin is chosen with equal probability. 3. the chosen coin is flipped and returned to the bag. The business goal is to predict the rate of heads in future trials (probably using Bayesian inference). This got me thinking - am I not over-complicating things? Isn't this process (observationally) equivalent to a single biased coin? Thanks in advance!


title: Species scaling methods in PCoA or CAP using vegan::capscale?

text: I want to study relations between sites categories and species abundances through PCoA or CAP using vegan::capscale. For doing so, I overlay species scores on my ordination. However, I am getting confused with the different scaling options and their interpretation. From what I understood: With scaling=1, arrow shows the direction from the origin for which sites have larger abundances for this species. With scaling=2, we want to analyze the correlations among species. Species that have a small angle between their arrows are expected to be strongly positively correlated. With scaling=-2, const = sqrt(nrow(dune)-1), we get correlations between species and axes. This comes from Jari Oksanen. I did compare the 3 different options (see codes below), differences seem to be only a matter of arrow length. Hence, by considering that the most important are the relative length of arrows (relative to each other), am I allowed to use scaling=-2 (species axes correlations) for both analyzing site- species associations and species-species correlations? One more questions, is there a well admitted threshold in the value of species-axes correlations, below which we consider that species are not correlated (I mean species don't differ in abundance across the sites, excluding the cases of non linear relationships). If yes, do this threshold change depending on the scaling method used and on either the ordination is constrained or not. Morover, if I want to overlay a vector for a quantitative environnemental variable, may I use also scaling=-2, from which correlation threshold?               library(vegan)     library(ggplot2)     library(grid)     data(dune)     data(dune.env)          dune=sqrt(dune)     mcap=capscale(dune~1,dist="bray") #PCoA          #sites scores     dims=c(1,2)     site=scores(mcap,display="wa",choices=dims)     site.env=cbind(site,dune.env)          #spider for management levels     dev.new();plot(mcap);coord_spider=with(dune.env,ordispider(mcap,Management,col="blue",label=F,choices=dims));dev.off()     coord_spider=as.data.frame(cbind(coord_spider[,],site.env))     names(coord_spider)[1:4]=c("MDS1","MDS2","MDS1end","MDS2end")          #species scores     #scaling1     cor.min=0.6 #below this threshold, arrows will be not plotted because      correlation is considered too much week     cor_sp=as.data.frame(scores(mcap, dis="sp", scaling=1,choices=dims))     cor_sp$cor=with(cor_sp,sqrt(MDS1^2+MDS2^2))         cor_sp$sup=FALSE;cor_sp$sup[cor_sp$cor>=cor.min]<-TRUE     cor_sp$lab=row.names(cor_sp)         cor_sp=cor_sp[cor_sp$sup==TRUE,]     cor_sp_s1=cor_sp          #scaling2     cor.min=0.6 #below this threshold, arrows will be not plotted because      correlation is considered too much week     cor_sp=as.data.frame(scores(mcap, dis="sp", scaling=2,choices=dims))     cor_sp$cor=with(cor_sp,sqrt(MDS1^2+MDS2^2))         cor_sp$sup=FALSE;cor_sp$sup[cor_sp$cor>=cor.min]<-TRUE     cor_sp$lab=row.names(cor_sp)         cor_sp=cor_sp[cor_sp$sup==TRUE,]     cor_sp_s2=cor_sp          #scaling -2, correlations between species and axes     #from Jari Oksanen at      https://stat.ethz.ch/pipermail/r-sig-ecology/2010-August/001448.html     cor.min=0.6 #below this threshold, arrows will be not plotted because      correlation is considered too much week     cor_sp=as.data.frame(scores(mcap, dis="sp", scaling=-2, const =      sqrt(nrow(dune)-1),choices=dims))     cor_sp$cor=with(cor_sp,sqrt(MDS1^2+MDS2^2))         cor_sp$sup=FALSE;cor_sp$sup[cor_sp$cor>=cor.min]<-TRUE     cor_sp$lab=row.names(cor_sp)         cor_sp=cor_sp[cor_sp$sup==TRUE,]     cor_sp_s3=cor_sp          #plot     mon.plot1=ggplot(data=site.env)+theme_bw()+        geom_point(aes(x=MDS1,y=MDS2,color=Management))# les sites          #add spider     mon.plot2=mon.plot1+     geom_segment(data=coord_spider,aes(x=MDS1,y=MDS2,xend=MDS1end,yend=MDS2end,colour=Management),lwd=0.5,alpha=1/3)+     geom_point(data=coord_spider,aes(x=MDS1,y=MDS2,colour=Management),cex=3,shape=19)          #add species scores as vector     #scaling1     mon.plot_s1=mon.plot2+ggtitle("Scaling 1")+     geom_point(aes(x=0,y=0),shape=21,fill="black",color="black",size=3)+#central      point        geom_segment(data=cor_sp_s1,aes(x=0,y=0,xend=MDS1,yend=MDS2),arrow =      arrow(length = unit(0.3,"cm")))+#arrows        geom_text(data = cor_sp_s1, aes(x = MDS1, y = MDS2, label = lab),      size = 3)#labels          #scaling2     mon.plot_s2=mon.plot2+ggtitle("Scaling 2")+     geom_point(aes(x=0,y=0),shape=21,fill="black",color="black",size=3)+#central      point        geom_segment(data=cor_sp_s2,aes(x=0,y=0,xend=MDS1,yend=MDS2),arrow =      arrow(length = unit(0.3,"cm")))+#arrows        geom_text(data = cor_sp_s2, aes(x = MDS1, y = MDS2, label = lab),      size = 3)#labels          #scaling -2, correlations between species and axes     mon.plot_s3=mon.plot2+ggtitle("Scaling -2")+     geom_point(aes(x=0,y=0),shape=21,fill="black",color="black",size=3)+#central      point        geom_segment(data=cor_sp_s3,aes(x=0,y=0,xend=MDS1,yend=MDS2),arrow =      arrow(length = unit(0.3,"cm")))+#arrows        geom_text(data = cor_sp_s3, aes(x = MDS1, y = MDS2, label = lab),      size = 3)#labels          #plot everything     grid.newpage()     pushViewport(viewport(layout = grid.layout(2, 2)))     print(mon.plot_s1,vp=viewport(layout.pos.row = 1, layout.pos.col = 1))     print(mon.plot_s2,vp=viewport(layout.pos.row = 2, layout.pos.col = 1))     print(mon.plot_s3,vp=viewport(layout.pos.row = 1, layout.pos.col = 2))


title: What does population actually refer to?

text: This is probably going to sound extremely stupid, but it's part of a larger question, about sample distribution of means, that I'm having problems understanding, so please bear with me. If I record the time between planes landing, and have 250 observations, are those 250 observations my "population", or does the term "population" refer to all the values, i.e. if I had sat there till the end of time recording? I'm assuming its the latter but just want to double check.


title: Sampling with or without replacement?

text: I don't know a lot about sampling methods. I have a large population of size 2,000,000. I used one of those sample size calculators. It says that I need sample size of approximately 10,000. I am trying to find the probability **p** of success for the population. It is not feasible for me to test all 2,000,000 members of the population. That is why I am sampling. I assume that the sample size calculator means a simple random sample without replacement. I have read that a simple random sample with replacement ensures that the covariance between two variables is 0 ~~, i.e., independent~~. When should one choose with replacement instead of without replacement? If we sample with replacement, then we are simply performing Bernoulli trials. I suppose this makes applying statistical tools (which?) easier. Again, sampling ignoramus here.


title: How to organise data and run nerlovian model

text: I am trying to estimate the price elasticity of supply for small scale farmers in Malawi and I have time series data for 34 years. I have two problems: First, the prices are very low such that when I take their logs am getting negative values, as a remedy i added 1 to each value and i got positive logs but i don't know how to interprete the result with the 1 that I added. Second i want to find out the specific way of estimating the model in `Stata`, I just used regress. I will appreciate your help


title: Football League coincidence

text: In the English Championship division, there are 24 teams, 8 of which have names starting with the letter B (e.g.Bolton) Tonight, (5th March 2013) all 24 teams in this division are playing each other. By a coincidence, the 8 teams starting with B are playing each other, i.e. 4 of the games involve these 8 teams (there are 2 teams per game!) What is the probability of this happening ?


title: PCA on linguistic variables: How far languages are?

text: I performed principal components analysis on continuous variables describing 16 languages. Using the first two axes, which explain 76% of variance, I need to calculate the distance between each pair of languages as appeared on the first two axes; to test in a Mantel test the correlation between distances in linguistic variables and geographic distances. Could anyone help me: how can I do that? cheers


title: Backfilling ARIMA data with exogenous variable

text: I have time series data for a set of cities that goes back for about 10 years. I also have the data at the state level for almost 30 years. There was an event that occurred about 20 years ago, that is captured in the longer, state level data, but not the city data, that I would like to investigate at the city level. What I think might be useful is to create some kind of ARIMA model that regresses the state data as an exogenous variable. If I were to do this, how do I use the model to backfill the city data such that it ends at the same point as where the actual city data starts? Is there already a canonical method of doing something like this? Thanks for any help you can give (literature references, R libraries, etc.)


title: What is a good resource that includes a comparison of the pros and cons of different classifiers?

text: What is the best out-of-the-box 2-class classifier? Yes, I guess that's the million dollar question, and yes, I'm aware of the no free lunch theorem, and I've also read the previous questions:   * Poll: What is the best out-of-the-box 2-class classifier for your application?   * and Worst classifier Still, I'm interested in reading more on the subject. What is a good source of information that includes a general comparison of the characteristics, advantage, and features of different classifiers?


title: Estimate confidence interval of mean by bootstrap t method or simply by bootstrap?

text: When estimating the confidence interval of mean, I think both the bootstrap t method and the nonparametric bootstrap method can apply, but the former requires a little more computation. I wonder what the advantages and disadvantages of bootstrap t over the normal nonparametric bootstrap are? Why? Are there some references for explaining this?


title: are there any statistically significant differences between my study population and a sub-group within it?

text: I have a total population of 31, of these I have a treatment finishing population of 20. I want to run a test on SPSS to identify if there are any significant differences between the 20 individuals that completed the treatment and the group as a whole. The variables I am trying to compare are age, weight, height and BMI. I have run descriptives for both the entire population and the sub group seperatley, I can see there is little difference between the groups but would like to know the statistical significance


title: Comparing the mean of several exponentially distributed variables

text: (Possibly related to Is there a nonparametric equivalent of Tukey HSD?) Given a set of several exponentially distributed variables (representing data which is modeled this way), I'd like to check whether one of them has a significantly higher mean than the others. If the variables were normally distributed, I would have used Tukey's test; however, this is not the case, and I'm looking for as simple alternative as possible. **Some background on the data** Alternative methods and markets from which we buy goods are compared every day. Data is aggregated on a daily basis (i.e. I only have prices and number of transactions per day, which translate to a mean). The distribution of funds for the next day is based on previous' day analysis: best method gets most of the funds, and the other methods share the rest. Typically, 5-10 alternative methods are compared, with 10,000 and more transactions per method per day. **Current Alternatives** Based on my limited knowledge, I came up with two possible ways of comparing methods:   1. Based on the referenced question above, a nonparametric versions of Tukey's test can be used.   2. Since only the highest-mean variable is compared to the rest of the variables, I think it can be compared to the second highest with _bonferroni correction_. i.e. use effective alpha = 0.05 / number of variables. This might be very conservative, but since the numbers are high, the results are usable (i.e. most days have significent results). **My questions**   1. Do the methods above make sense as solutions to the presented case? specifically, method (2) is simple and yields usable results, but I'm not sure it's statistically sound.   2. I'm interested in any other suggestions for performing this analysis. As noted, I'm more interested in simplicity and generalization, and less on getting the best significance values.


title: Distribution help

text: In the following problem, check that it is appropriate to use the normal approximation to the binomial. Then use the normal distribution to estimate the requested probabilities. It is known that 85% of all new products introduced in grocery stores fail (are taken off the market) within 2 years. If a grocery store chain introduces 68 new products, find the following probabilities. (Round your answers to four decimal places.) (a) within 2 years 47 or more fail (b) within 2 years 58 or fewer fail . (c) within 2 years 15 or more succeed . (d) within 2 years fewer than 10 succeed


title: E-step in EM-algorithm using MAP estimate (mixed Markov models), what does it calculate?

text: I am trying to grasp what exactly is "estimated" in the E-step of the algorithm. According to all definitions, in E-step the "conditional expectation values , or posterior probabilities of the hidden variables" are computed, using the Bayes formula ( posterior probability=prior( or marginal prob.) x likelihood/probability of evidence). Now, my question is, are these posterior probabilities, that I "estimate" or "compute", actually numbers or are they functions with respect to the current iteration's parameters ( and they are the ones I need to plug into M-step)? Where the probability of evidence would be likelihood, calculated with the previous argmax , calculated in previous M-step.


title: 2 WAY ANOVA if groups are missing

text: I like to set up a 2 way repeated measures anova. I compare 4 models with 4 ages and see if the results differ, this should result in around 16 groups, since 4*4=16 (so far so good) However from some groups I don't have the data so I still have 2 factors to test for but not longer 16 groups but 15-14. This is not due to mistakes but due to the fact that the groups under this condition are just bad and don't meet standard recommondations. If however I'll try to do 2 way repeated measurments including my missing groups, I can't execute it. Is there a way to handle this missing groups? Greetings


title: Deep dive into recall, precision and threshold calculations in R

text: I have a textual classification problem that consists of two categories- zero an one. Up until now I tried solving it by creating a Document Term Matrix, and to run it through SVM (using RTextTools package). Here's a code snippet: (R)               models <- train_models(container, algorithms=c("SVM"))     results <- classify_models(container, models)     analytics <- create_analytics(container, results)     View(summary(analytics))          >>ALGORITHM PERFORMANCE          >>SVM_PRECISION    SVM_RECALL    SVM_FSCORE      >>         0.64          0.63          0.63       My questions are as follows:   1. Why are all the predicted values in the result matrix between 0.5-1? isn't it supposed to be 0-1?    2. How can I see (in R) under which **_threshold_** are these precision and recall values being calculated? How can I change this threshold to get different values?   3. How exactly are the recall/precision scores being calculated? I mean, supposed we have _theta_ to separate that all scores above it are of class 1, and the rest are 0. Are these output scores just a simple mean of the two classes precision/recall? Is it a weighted average?   4. How can I create in R two different thresholds values for each class (with what's left in between labeled as "unidentified")?


title: Fitting random intercept and slope in lmer for lme4.

text: If I would like to fit random intercept and slope and if I write it as (color|writer) compared to (1+color|writer), are they the same?


title: Filter out "abnormally large" values from a list of data

text: For example, the lists can be something like: $$\{1.123213, 5.154543, 2.134121, 7.34534, 12.223432, 8.16571, 100.45645, 222.423\}$$ I want to remove $\{100.45645, 222.423\}$. $$\{232.123213, 323.154543, 232.134121, 123.34534, 222.223432, 8545.16571,\\ 4335.45645, 1222.423\}$$ I want to remove $\{8545.16571, 4335.45645, 1222.423\}$. The size is not fixed, and the range may vary. Are there any ways to remove these "abnormally large" values from my lists?


title: How to get a reciprocal plot from x and p(x)?

text: Let's say that I have a picture like the one in figure. This look like an histogram, but actually is not, in the sense that I did not produce the figure starting from a dataset of values. In fact, I only have the probability (the y axis) for each bar. I also have the center point of each bin, and the bin width of each bin. Notice that the graph is normalized such that the area under the bins is equal to one. ![enter image description here](http://i.stack.imgur.com/CQWGd.jpg) What I want to do now, is to calculate the reciprocal of this probabilities. Usually I would use this formula: $g(y) = \frac{ 1 }{ y^2 } f\left( \frac{ 1 }{ y } \right) $ But I cannot do that now, since I do not want to assume any f(x). Having an histogram, I would just plot the (1./dataset) to get the reciprocal (in histogram form), and then I would normalize (assuming that the original dataset was not already normalized like the one in the example), so to have a good approximation of the rate distribution. However, in this case, I don't have the dataset. I only have the y points (the p(x), actually) and the x (bin center for each bin). I have been adviced to use the following formula (in MATLAB): plot(1./x, y) And this will give me a function with the same shape as the reciprocal function, but not normalized. I _think_ I could take this result and normalize it, and that should work fine. However, I wonder if there is another method that uses only the x and the p(x) in order to get the reciprocal values. Thank you.


title: EnergyThe meaning of these term in Image Processing?

text: Im newbie in Image Processing field, and I got many term that wierd. There are :   1. Moments   2. Local and global characteristic   3. Energy of Image   4. Color Temperature What each of them represent in Image??


title: How to solve this example using Poisson Distribution?

text: Particles are suspended in a liquid medium at a concentration of 4 particles per mL. A large amount of the suspension is thoroughly agitated, and then 3 mL are withdrawn . Let $X$ be the number of the particles in the 3 mL. Answer the following: a- The distribution of X is a Poisson distribution, what is its parameter ($\lambda$)? b- Find the probability that no particle is withdrawn from the 3 mL. c- Find the probability that one particle is withdrawn from the 3 mL.


title: Why is the probability zero for any given value of a normal distribution?

text: I noticed that in the Normal distribution, the probability $P(x=c)$ equals zero, while for the Poisson distribution, it will not equal zero when $c$ is a non-negative integer. My question is: Does the probability of any constant in the normal distribution equal zero because it represents the area under any curve? Or it is just only a rule to memorize?


title: Approach to solving for the probability that a photo contains an object based on tags

text: I've studied machine learning and statistics only very briefly. I've used linear regression to solve problems with fixed sets of variables, but I'm not sure how to approach the following problem. Given a bunch of tagged photos, I'd like to use the tags to estimate the probability that a specific object is contained within the photo. I have access to a wealth of data which has already been classified as either containing the objects we're looking for or not (however note that this is not a classification problem; the probabilities are being used to rank photos). For example. If a photo is tagged with _coke_ and _drink_ , there is a high probability it contains coca-cola. However if a photo is tagged with _coke_ and _crack_ , there is a high probability it does not contain coca-cola. If a photo is tagged with either _crack_ or _drink_ separately, that doesn't tell us much about what might be in the photo. How would one go about building a hypothesis formula for this problem?


title: Interpreting coefficients in a logistic regression model with a categorical variable having more than 2 levels

text: There is quite some content online interpreting odds in a logistic model with a dichotomous predictor. My problem is understanding coefficients when there are more than 2 levels for a categorical variable. How do you _define the odds_ then?               Data:      X is a single categorical predictor with 4 levels: teenager, adult, mature, senior.            Y: 1=smoking, 0=non smoking.          LR:  We use n-1 dummy variables.            I chose adult as the reference bin as it had the highest concentration. (ok??)          ________ | Intercepts   | p     adult    |    -4.3801   | 0     teenager |    -0.32456  | 0     mature   |     1.45119  | 0     old      |    -0.9891   | 0      **Interpreting the coefficients** Teenager: Teen is less likely to smoke (w.r.t adult?). In fact, a teen is 28% (exp-0.32456 -1) less likely to smoke THAN AN ADULT. Is odds of teenager smoking always mentioned against the reference group? Mature: Matures is more to smoke (w.r.t adult?). In fact, a mature is 326% more likely to smoke THAN AN ADULT. Is odds of mature smoking always mentioned against the reference group?


title: Probability of event happening when data is aggregated with many independent events over the course of time

text: Let's say you have $X$ coins, each with a differing probability of landing heads (e.g. coin 1 has 10% chance of landing heads, coin 2 has 20% chance of landing heads, etc.). Now, let's say that you flip coin $i$ coin $Y_i$ times (each coin has a differing amount of flips). We know how many times each coin was flipped. Now let's say you do this every day for a really long time and record that info. Example, with 2 coins, on day 1: 50 flips total, 30 heads total, coin1 was flipped 20 times and coin2 flipped 30 times. Day 2: 80 flips, 66 heads total, coin1 60 flips, coin2 20 flips. What we know: total flips, total heads, how many times each coin was flipped. Given that, is there a way to determine an approximate probability that the given coin will flip heads? In the above example, coin1 has a 100% probability of heads and coin2 has a 33% probability.


title: Positive coefficient but negative marginal effect in mlogit

text: Is it plausible to have a positive coefficient with a negative marginal / impact effect after running multinomial logit model?


title: Temporal trends with missing values

text: I need to calculate the temporal trends for some climate variables with missing values. For example, last frost days defines as the last day of year with minimum temperature less than 0C. However, there are no any frost days in some years. My data look like:               lfd <- c(NA, NA, NA, NA, NA, 190, 192, 189, 200, 185, 205, 203, 200, 207, NA, NA, 205)     years <- seq(1957, length.out = length(lfd))      Now I use the linear regression in R (function lm) to calculate temporal trend. It seems the results are unreasonable for datasets with many missing values. How could I calculate the temporal trends with missing values? Thanks for any suggestions.


title: Designing genomic search strategies for non-coding RNAs using BLAST

text: I am interested in developing a computational search of some types of non- coding RNAs (ncRNAs) on a new genome. I have some questions related to the nature of my experiment, because it is an exploratory study, and I can generate different results by applying different BLAST search strategies. The final output consists of different counts of each type of ncRNA for each BLAST strategy. How I can design good search strategies, taking into account that it is the first search on this new genome? References will be appreciated. Thanks!


title: How to determine if a risk indicator is really a risk indicator

text: Let's say I ask a robot a question: "Why is partial differential equations no fun at all?" Suppose the robot is stupid and categorizes my question inside the class of "cats" so it replies with: "Cats are fluffy and awesome." So I reword my question over and over and sadly the robot keeps hitting the "cats" class. Sequential hits is an example of a risk factor. It indicates that perhaps something is wrong with the "cats" class in the robot's brain. Maybe the word "PDEs" appears as sample text under the class to draw from. Suppose I have a collection of classes (e.g. "cats", "dogs") that all share the same risk factors(e.g. "sequential hits","bad user scores",etc). Each risk factor has its own weight which I can adjust from 0-10 where 0 essentially means that I can completely drop the risk factor and 10 is the most important. For each class, I can multiply the occurrence of each risk factor by its respective weight and sum them all to obtain a risk score for the class. Thus, each class has a risk score and I can focus on fixing the classes with the highest risk scores. My question is this: How can I optimally choose the weights of the risk factors to give me what are truly the riskiest classes? This is a pretty open question. I am looking for ideas, papers, books, or statistical techniques that might be helpful to this endeavor. Much thanks.


title: Confusion related to posterior gaussian process

text: I was reading this `[paper][1]` related to sparse online gaussian processes. However, I didn't get how the denominator in the equation 1 was derived?![enter image description here](http://i.stack.imgur.com/vr3dG.png) It was supposed to be P(D). I didn't get how the denominator was obtained. Further, the paper is saying that the posterior may not be gaussian. I don't know why it is so?


title: Analysis of ordinal data with multiple assessments per unit... the right approach?

text: Suppose I have a dataset with the following information:   * N objects, each of which can be rated   * The number of ratings per object is variable (raters are unidentified), e.g., there are 30 objects each with a unique id   * Each row has the predicted rating, the actual rating, and the difference between the two, which is an ordinal value, say between 0 and 10 How do I measure the consistency of ratings per object and for the whole dataset? Basically, there is a distribution of "differences" per object. What are the various ways to approach this problem? Thank you.


title: Formula suggestion for a scoring system, weighted solution

text: I'm trying to validate a series of words that are provided by users. I'm trying to come up with a scoring system that will determine the likelihood that the series of words are indeed valid words. Assume the following input:               xxx yyy zzz      The first thing I do is check each word individually against a database of words that I have. So, let's say that `xxx` was in the database, so we are 100% sure it's a valid word. Then let's say that `yyy` doesn't exist in the database, but a possible variation of its spelling exist (say `yyyy`). We don't give `yyy` a score of 100%, but maybe something lower (let's say 90%). Then `zzz` just doesn't exist at all in the database. So, `zzz` gets a score of 0%. So we have something like this:               xxx = 100%     yyy = 90%     zzz = 0%      Assume further that the users are either going to either:   1. Provide a list of all valid words (most likely)   2. Provide a list of all invalid words (likely)   3. Provide a list of a mix of valid and invalid words (not likely) As a whole, what is a good scoring system to determine a confidence score that `xxx yyy zzz` is a series of valid words? I'm not looking for anything too complex, but getting the average of the scores doesn't seem right. If some words in the list of words are valid, I think it increases the likelihood that the word not found in the database is an actual word also (it's just a limitation of the database that it doesn't contain that particular word). NOTE: The input will generally be a minimum of 2 words (and mostly 2 words), but can be 3, 4, 5 (and maybe even more in some rare cases). **EXAMPLE 1:** Say the following scores:               xxx = 100%     yyy = 100%     zzz = 0%      The average is 66.66%. But since already two words in the list exist in the database, the chances are that `zzz` is also a real word. Scoring this series of words as just 66.66% seems low.


title: Bayesian complexity tests

text: I have seen an apparently interesting quantity called the Bayesian complexity, which in this case is defined as $-2(\overline{\ln\mathfrak{L}} - \ln\mathfrak{L}_\mathrm{max})$. I am not sure what this is used for. I.e., what sorts of tests can be made on this quantity and what does it imply?


title: Outlier Removal in Non-Normal Data

text: I am currently working on a project in which I have to eliminate outliers from non-normally distributed data sets. The data sets are subsets of a fairly large database (order of millions of observations) segmented into groups ranging from 30-200+ observations based on the type of product being looked at. The data is pricing data which could be abnormal in the data set due to the possibility of rush orders and other types of abnormal circumstances under which purchases are made which can drive prices above or below a typical range). The difficulty in trying to identify these outlier values is that not every group of observations follows a normal distribution so tests for outliers based on normality are not accurate in testing for outliers in this type of data. My current approach to this problem is to measure the standard deviation from the mean to flag all values less than three deviations from the mean of the group, this is used to flag the observations for administrative review later. then I measure skewness using the Pearson's coefficient of skewness. If the value of this statistic is beyond a certain threshold ±.4 then I get rid of the 1st percentile (if skew is lower than -.4) or I remove the 99th percentile of data (if the skew is greater than .4). This measurement of skew is repeated with a new measurement of skew each iteration with the same type of observation removal based on the parameter listed in the previous paragraph. The removed values are then given over to administrative review to determine whether they are true outliers with special attention being paid to any values within 3 standard deviations. is this a good approach to removal of outliers in this case? I don't have much experience with outlier detection in data not belonging to a normal distribution. so the steps listed explicitly are   1. Calculate deviation from mean for each observation   2. Calculate skewness of the pool of observation if skew>.4 remove 99th percentile observations if skew<.4 remove 1st percentile observations   3. values removed are then put under review to determine outlier status with special attention paid to values falling within 3 standard deviations of mean.


title: When should I use each of these methods to calculate correlation?

text: I am using R for data analysis. R provides a `corr` function for calculating the correlation. This function provides three different approaches/algorithms to estimating the `corr` which are Pearson, Spearman and Kendall. When should I use each of each of these methods? What factors determine which method should be used?


title: How to deal with time series with revisions?

text: I manipulate in R time series that can present **revisions**. I use in general the package 'xts' but it doesn't support revisions. Is there a package that deals with that? Alternatively, is there a best practice to store and manipulate revisions?


title: How to calculate mean on group statistics

text: ![enter image description here](http://i.stack.imgur.com/qTI9v.jpg) How can I calculate the mean of the first schema? I don't have the individual values of the sample. Std. Deviation (B) is 0.6441.


title: How to reliably estimate a random sum of unreliable elements?

text: This is an ill posed problem I'm faced with at work. I have observations $X_i$, $i = 1,2,3,...N > M$. I need want to compute $S = \sum_{i=1}^M X_i$. Unfortunately, $M$ is random, and extremely noisy. It turns out, I have some rough "goodness" measures, $w_i$ for each of the $X_i$. That is, $w_i$ low, means I don't trust $X_i$ very much, $w_i$ high, means I trust $X_i$ a lot. Moreover, "generally" as $i$ gets large, the $X_i$ get small. (There's a lot of dependence between $M$ and the $X_i$) In this application, it turns out that the _average_ of the $X_i$ can very reliably (i.e. with relatively small standard deviation) be computed by: $A = \frac{\sum_{i=1}^M w_i X_i}{\sum_{i=1}^M w_i}$. We know the standard deviation of $A$ is relatively small by running a number of independent trials, and computing the sample standard deviation. But, I similarly need a "good" estimate (i.e. small standard deviation) of the straight sum of the $X_i$ and the obvious approach: $M * A$ has far too much variance for our application. Understand this is pretty vague, but any thoughts on how to proceed?


title: Which correlation analysis is appropriate for binary and ordinal variables?

text: Violin players were asked to rate 5 violins from least to most rich on a `[0 1]`, `0.05-increment` scale, `most always = 1` and `least always = 0`. They were then asked to set a "limit of acceptability," i.e., to set a point along the scale above which violins are considered acceptable in terms of richness. By setting `acceptable = 1` and `not acceptable = 0`, a binary variable is obtained for each player, in addition to their rating. So, as an example, for a player I have the variables `rating = [0.3 0.45 0 1 0.85]` and for a `limit = 0.45`, `acc = [0 1 0 1 1]`. I use Lin's concordance correlation coefficient to explore consistency within and between individuals for the case of ratings. However, this is not possible for the "acceptable" variable: when a player sets their `limit = 0`, then `acc = [1 1 1 1 1]` and hence the presence of zero variance becomes an issue. I am trying to figure out a good, interesting way to analyze the binary data, perhaps compare consistency in ratings with consistency in acceptability limits, or more in-depth relations between the two variables. I am not very experienced in statistics, so any suggestions and ideas are welcome. I primarily work with MATLAB.


title: How would you preprocess data for SVD?

text: I am computing SVD on a matrix which is the empirical version of $E[XY^{\top}]$ for some $X \in \mathbb{R}^{m \times 1}$ and $Y \in \mathbb{R}^{n \times 1}$. I am wondering if there are standard ways to preprocess $x_1,\ldots,x_l$ and $y_1,\ldots,y_l$ before doing that (other than subtracting the mean and dividing by standard deviation). This is related to the question here: "Normalizing" variables for SVD / PCA.


title: Thinking about how adding extra variables affect linear regressions

text: I have two models. First, $$z=\alpha_1x+\alpha_2+\epsilon$$ and second, $$z=\beta_1x+\beta_2y+\beta_3+\epsilon.$$ What factors affect whether $\alpha_1\approx\beta_1$? Note that I don't care how accurately either model predicts $z$. I only care about how similar the coefficient of the $x$ term is in both models. Assume both $x$ and $y$ actually predict $z$ and don't interact.


title: Maximum likelihood estimation when parameters are functions of another data series

text: We have two time series: $X_t$ and $R_t$, and a model saying that $R_{t+1} = (\mu(X_t) - \frac{1}{2}\sigma^2(X_t))\Delta T + \sigma(X_t) \sqrt{\Delta T} \epsilon_t$, where $\Delta T$ is given constant and $\epsilon_t$-s are independent normally distributed with zero mean and unit variance. Further we assume that the functions $\mu(x)$ and $\sigma(x)$ are linear for simplicity. I would like to use some standard method (MLE comes to my mind) to estimate parameters of functions $\mu(x)$ and $\sigma(x)$, but I am not sure how to do this. I would be grateful for detailed answers, because I am not really experienced with statistics.


title: How to interpret statistics (Emax, D, U, Q, B, etc.) of bootstrap validation of logistic regression

text: I'm only a linguist, so my knowledge of statistics is very basic. I fitted a logistic regression model with R (with `lrm(formula, y=T, x=T)`), and when I use the option `validate(lrm)`, I get some statistics I don't really understand.               index.orig   training   test optimism index.corrected     n     Dxy           0.5984   0.6112  0.5461   0.0651          0.5333 40     R2            0.3258   0.3676  0.2929   0.0747          0.2511 40     Intercept     0.0000   0.0000 -0.0105   0.0105         -0.0105 40     Slope         1.0000   1.0000  0.8427   0.1573          0.8427 40     Emax          0.0000   0.0000  0.0399   0.0399          0.0399 40     D             0.2713   0.3176  0.2394   0.0782          0.1931 40     U            -0.0177  -0.0177  0.0092  -0.0269          0.0092 40     Q             0.2890   0.3353  0.2302   0.1051          0.1839 40     B             0.1864   0.1772  0.1972  -0.0201          0.2064 40     g             1.4632   1.6642  1.3460   0.3182          1.1449 40     gp            0.2840   0.3011  0.2703   0.0308          0.2532 40      I don't really understand most of that. I think `R2` and `Dxy` are supposed to be statistics of how good the predictors are, but I'm not sure how I should interpret the values, does the corrected `Dyx = 0.651` mean that there is a strong correlation, while the corrected `R2 = 0.0747` means that the correlation is very weak? I think the model is overfitted, but I'm not sure if I'm right. Also, the other statistics are totally strange to me. What are `Emax, D, U, Q, B, g`, and `gp`?


title: Implementation of Gaussian process

text: i try to implement a simple gaussian process regression in java. I almost got every step from the book http://www.GaussianProcess.org/gpml . With my implementation of algorithm 2.1 on page 19 i'm able to produce results but: I have the strange behavior that all predicted datapoints where i do not have any target tend to be of zero value. Does anyone know this behavior or can explain where i could have taken a mistake? In my plot the red line marks the predicted values, black are the targets. Here is the prediction and covariance java code. The parameter array includes length scale and noise, I work with the idexes of the arrays as x inputs for the covaiance function:               private double generateCOV_RADIAL(int i, int j) {              double covar = Math.pow(parameter[1], 2.0) * Math.exp(-1.0 * (Math.pow(i - j, 2.0)/(2 * Math.pow(parameter[0], 2.0))));              if(i==j)                 covar += Math.pow(parameter[2], 2);               return covar;     }     public void predict(int predictFuture) throws IllegalArgumentException, IllegalAccessException, InvocationTargetException{              RealMatrix y = MatrixUtils.createColumnRealMatrix(this.targets);              double mean = StatUtils.mean(this.targets);              for(int i = 0; i < this.targets.length;i++)         {             targets[i] -= mean;         }                   RealMatrix K =  MatrixUtils.createRealMatrix(cov);                   //identity matrix for I         RealMatrix k_eye = MatrixUtils.createRealIdentityMatrix(cov.length);                   //choleski(K + sigman^2*I)         CholeskyDecomposition L = null;         try {             L = new CholeskyDecompositionImpl(                     K.add(                             k_eye.scalarMultiply(Math.pow(parameter[2], 2))                             )                     );         } catch (NonSquareMatrixException e) {             e.printStackTrace();         } catch (NotSymmetricMatrixException e) {             e.printStackTrace();         } catch (NotPositiveDefiniteMatrixException e) {             e.printStackTrace();         }              //inverse of Ltranspose for left devision         RealMatrix L_transpose_1 = new LUDecompositionImpl(L.getLT()).getSolver().getInverse();         //inverse of Ltranspose for left devision         RealMatrix L_1 = new LUDecompositionImpl(L.getL()).getSolver().getInverse();                        //alpha = L'\(L\y)         RealMatrix alpha = L_transpose_1.multiply(L_1).multiply(y);                   double L_diag = 0.0;              for(int i = 0; i < L.getL().getColumnDimension();i++)         {             L_diag += Math.log(L.getL().getEntry(i, i));          }              double logpyX = - y.transpose().multiply(alpha).scalarMultiply(0.5).getData()[0][0]                         - L_diag                         - predictFuture * Math.log(2 * Math.PI) * 0.5;                   double[] fstar = new double[targets.length + predictFuture];         double[] V = new double[targets.length + predictFuture];              for(int i = 0;i < targets.length + predictFuture;i++)         {                  double[] kstar = new double[targets.length];                  for(int j = 0; j < targets.length;j++)             {                 double covar = (Double)covMethod.invoke(this,j,i);                 kstar[j] = covar;             }                  //f*=k_*^T * alpha             fstar[i] = MatrixUtils.createColumnRealMatrix(kstar).transpose().multiply(alpha).getData()[0][0];             fstar[i] += mean;             //v = L\k_*             RealMatrix v = L_1.multiply(MatrixUtils.createColumnRealMatrix(kstar));                  //V[fstar] = k(x_*,x_*) - v^T*v             double covar = (Double)covMethod.invoke(this,i,i);                  V[i] = covar - v.transpose().multiply(v).getData()[0][0] + Math.pow(parameter[2],2);         }              this.predicted_mean = fstar;         this.predicted_variance = V;     }      Thank you ![GP Example](http://i.stack.imgur.com/hfjBL.jpg)


title: Need help in creating an equation of meeting quality

text: I want to create an equation that consists of two independent variables "duration" and "data_transferred", and one dependent variable "meeting_quality" This is what I want: meeting_quality would be good, if more data is transferred in less duration. meeting_quality would be bad, if less data is transferred in more duration. What should be the equation presentation of above sentence? Thanks in advance. In fact, I want to use the values above equation in regression analysis. I would really be grateful for any further advice.


title: How to present rows which can have several categories?

text: I have a spreadsheet where the content is organized like that:               content  | category1 | category2 | category3 | ...     _________|___________|___________|___________| ...     content 1|     x     |           |     x     | ...     content 2|           |     x     |           | ...     ...      How can I present this information more effectively? For example, is there a Javascript library I could use or something like that, where we could choose to show only the rows which are associated to some categories, or set a category to a row, etc? Thanks in advance


title: Multiline names (x labels) in boxplot

text: In `R`, I would like to give multiline names to some data sets in `boxplot`, like this:               boxplot(rnorm(10), rnorm(10, mean=2), names=c("Normal", "Shifted\n*"))      Here, the names seem to be aligned with their bottom lines, which causes the text to overlap with the axis. How can I have multiline names which simply extend to the bottom instead of upwards?


title: Bubble charts with pie charts as bubble

text: I wanted to create a bubble chart which means a scatter plot with a separate variable indicating diameter of bubble. Is there any way to use bubble then as pie chart? I would really like to do this in MS Excel but I just don't know whether it is possible. * * * I have another problem with the code proposed earlier in the post. Now that the first problem appears to be solved there is a problem in the line `Set rngRow = Range(ThisWorkbook.Names(“PieChartValues”).RefersTo)` and I tried to replace PieChartValues with the range from the first row to the last in which the pie charts values ( the values from which the pie chart is generated) - like "A1:G12" and all the data are in that range. The tab in the sheet I am working in is called "B" without quotes should I write `Set rngRow = Range(ThisWorkbook.Names(“B”).RefersTo)` * * * i tried to implement the code of a link proposed earlier in VBA this is               Sub PieMarkers()          Dim chtMarker As Chart     Dim chtMain As Chart     Dim intPoint As Integer     Dim rngRow As Range     Dim lngPointIndex As Long          Application.ScreenUpdating = False     Set chtMarker = ActiveSheet.ChartObjects(“chtMarker”).Chart     Set chtMain = ActiveSheet.ChartObjects(“chtMain”).Chart          Set chtMain = ActiveSheet.ChartObjects(“chtMain”).Chart     Set rngRow = Range(ThisWorkbook.Names(“PieChartValues”).RefersTo)          For Each rngRow In Range(“PieChartValues”).Rows          chtMarker.SeriesCollection(1).Values = rngRow     chtMarker.Parent.CopyPicture xlScreen, xlPicture     lngPointIndex = lngPointIndex + 1     chtMain.SeriesCollection(1).Points(lngPointIndex).Paste          Next          lngPointIndex = 0          Application.ScreenUpdating = True          End Sub      when I execute the code I get an error pointing to the line Set chtMarker = ActiveSheet.ChartObjects(“chtMarker”).Chart does anybody have an idea why? I have an object a pie chart named exactly chtMarker so the object should be there


title: Difference between two i.i.d Laplace distributions?

text: What is the PDF of the difference of two i.i.d Laplace distributed random variables? I know that the difference of two i.i.d Normal variables is still the Normal distribution. Since the properties of the Laplace distribution are similar to the Normal distribution, I am guessing that the difference is also the Laplace distribution. I have tried to solve this using Mathematica. Neither `PDF` nor `Integrate` functions give a correct answer. Thank you!


title: Comparing two measurement devices without a criterion measurement

text: I am comparing energy expenditure outputs from two devices that track physical activity in a free-living setting over 10 days. Doubly labeled water is the accepted criterion measurement for this, but 1) This technique only gives total energy expenditure and 2) I did not use this in my study... Both of the devices I used have previously been validated against DLW. These devices give a lot more information than DLW, namely time spent at specific PA thresholds and number of 'bouts' at each threshold. I would like to assess the level of agreement between the two devices. I know that Bland-Altman plots are commonly used to assess agreement between measuring devices but I am struggling to decide which test to use without a reference/criterion measurement. I would guess that a basic correlation would provide some information but was wondering if anybody could offer some more information and/or guidance for a more sophisticated statistical anaylsis. Maybe a regression / sum of squares would be appropriated? Any help would be much appreciated.


title: Chi Square Analysis Percentages

text: I've used chi-square analysis in school, but was absent for the lesson, and I never quite got the hang of it. I've been successful with it in the past, but now I'm trying to use percentage values, and I'm not getting the answers I'm expecting. I'm doing a survey and comparing my findings with national averages, but when I use percentages (such as the percent of participants between the ages of 21 and 29) I get understandably tiny numbers (2.3% is 0.023). If I'm using percentages for the expected value, should it matter? Ex: Chi² for 21–29 age group is $(.084 - 0.1642)^2 / 0.1642 = .039187$ The total chi-square value is 0.1368, which is much smaller than the 0.5% chi squared table value for 8 degrees of freedom, 17.535, but the numbers seem way off. Do I have to find the whole numbers?


title: Identify sub-dimension correlation between 2 variables

text: If we have two variables A={20 values} and B={20 values} and we want to measure the correlation between these two variables. Lets assume that the first n values are highly correlated but the remaining m values are not correlated. The overall correlation between A and B may not represent the importance of the relationship. Imagine that there is only weak correlation in some values, but several values show really good correspondence between A and B. Traditional correlation coefficients will not capture such relationship. Is there any measure or statistical method that can detect such relationship between A and B?


title: SPM for spatial analysis

text: Do you know of any study where researchers have used Statistical Parametric Mapping for **spatial data only**? I found a lot of studies that used it for time-series analysis, but does this analysis work for data that don't share a time dimension?


title: Practical interpretation for $u_t = \log(x_t) - \log(x_{t-1})$

text: I have a time series $x_t$. If I use the transformation $u_t = log(x_t) - log(x_{t-1})$, my new time series $u_t$ has properties of white noise (random). I wonder whether there is any practical interpretation for $u_t$? Thanks for your help.


title: Interpreting numerical value of standard error of the mean

text: I am unclear about how to interpret the value of the Standard Error of the Mean (SEM) directly. For example, when a mean is reported as 5.00 + 0.50SEM, how do you directly relate the 0.50 to 5.00? To quickly run through the basic theory concerning the standard error:   1. The standard deviation (SD) is a measure of dispersion around the mean    2. The SEM is the SD of the sampling distribution for the sample mean   3. This sampling distribution is derived from the means of an infinite number of samples from a statistical population and is normally distributed according to the Central Limit Theorem   4. In a normal distribution, 68.3% of (randomly selected) values fall within +1SD, 95.4% +2SD and 99.7% +3SD   5. The SEM decreases with increasing sample size QUESTION 1: Given points 2) and 4) is it correct to interpret the SEM in the same way as SD as a descriptive statistic? That is, for a sample with mean 5.00 and SEM 0.50, is it correct to conclude the true population mean lies between 4.50 and 5.50 with probability 68.3%? QUESTION 2: For a given statistical population being sampled, does the sampling distribution change with sample size (given point 5)? i.e. should point 3) be: "...derived from the means of an infinite number of samples of a given size from a statistical population..."? QUESTION 3: Since the SEM is not calculated directly but estimated from the SD of a sample, what effect does departure from a normal distribution of the sample have on calculation of the SEM? Put another way, if a sample has for example a highly skewed distribution, calculating the SD is inappropriate (it is a property of a normal distribution only), so can the SEM still be estimated reliably from the SD?


title: Possible outcomes of approximate profile-likelihood estimator (APLE) for spatial autocorrelation

text: I've been working with spatial autocorrelation for a while and now I'm trying to move from more traditional estimators such as Moran's I or Geary's C to the new APLE estimator. I read Li's papers on APLE and also the reference on R-project. But I am still missing something very basic, as I still didn't understand something: the results of the APLE statistic are said to have a closed form but I didn't figure which one it is. For instance, while Moran's I can gom from -1 to 1, I couldn't figure which are the possible boundaries of APLE statistics. More interesting, I've been running tests in my older data, using the spdep package in R, and it seems to me that APLE can go up to 2. Any one has any hint on that? Please, any insights here would be very much appreciated.


title: Resampling with replacement permutation analysis

text: My question might be easy for most of you to answer. I am starting to learn statistics and coding in R so my questions are on the basic level. I have a dataset of (two groups, replicates). My data is split based on the groups so I have 24 samples in group 1 and 20 samples in group 2. My data has replicates. So each set has 4 replicates, hence I have 6 sets in group 1 and 5 sets in group 2. Hence I have assigned indices to them to make it easier during permutation (indices from 1-11). What I want to do now is a routine permutation analysis to obtain the test statistic. I am using non paramteric method with resampling with replacement. My doubt\problem in R coding is that I have to pool the data together and then resample the groups. When I try to do this, **I have to make sure I maintain the sample size for respective groups (that is after resampling the group lables, my new dataset should still contain 6 sets (24 samples) in group 1 and 5 sets (20 samples) in group 2.** I am unable to achieve the latter. I tried something like this in my code:               n_resamples=100          permdat=function(dat){     for (i in 6:nc)     {      for (j in 1:n_resamples)        {             inds=unique(dat$ind)          cc=dat$grp[dat$ind==inds]         pcc=sample(cc,replace=T) ##for resampling with replacement         #m<-list()         pdat=dat          for(i in 1:length(inds))     {     pdat$grp[pdat$ind==inds[i]]=pcc[i]     #browser()     }      This does not produce the desired result. Kindly help with R codes. Thanking you all for your help. :) Regards Ap


title: Observed vs predicted values from a logit model

text: I have a logit model and am trying to understand and compare the predicted and observed values generated by the model. Let's say data set had 100 values and I generate all the predicted probabilities, and then I find the actual probabilities from the data set. If I'm comparing the predicted vs observed values, I'm thinking there are two ways to do it. One is to do it value by value, while the second would be to group by the 'predicted probabilities.'               Method 1:          x_value  pred_val   obs_val     100       0.30       0.34     102       0.33       0.36     104       0.35       0.37     106       0.40       0.40     ...      I'm also thinking there has to some way to aggregate these values. So I'm thinking of aggregating all x values where the predicted probabilities is between 10 to 20% percent, then find the avg predicted value from that range, followed by the predicted value for that range.               Method 2:          Pred_probs       pred_val    obs_val     10 to 20% vals     0.10        0.11     21 to 30% vals     0.12        0.16     31 to 50% vals     0.15        0.30      What I'm wondering is:   1. When there are a large number of data points, what use is having a list of the predicted and observed values for any given value of x?   2. Does it ever make sense to do something as identified in 'Method 2'?


title: Easy way to test usefulness of data for stock market analysis?

text: I have a lot of data (gigs) that may be useful in predicting equity prices. I can import these as a series of features (columns) in a table where the companies are rows. I have time series information too. I have some machine learning experience but no experience as a trader. Is there some software or platform where I could easily import my data and it could backtest/forward test my data to see if it's useful? I understand that any machine learning system that's out there in public won't outperform the market enough to cover brokerage fees, but given my data, there's a chance it would. So I probably don't need state of the art in machine learning but I'd like to find a solution where I don't spend the next 6 months learning the stock market. I'd rather spend that time getting feedback on, and iterating on the input data, because as a data-oriented developer, that's where I can add value. Any help, much appreciated.


title: Why doesn't k-means give the global minimum?

text: I read that the k-means algorithm only converges to a local minimum and not to a global minimum. Why is this? I can logically think of how initialization could affect the final clustering and there is a possibility of sub-optimum clustering, but I did not find anything that will mathematically prove that. Also, why is k-means an iterative process? Can't we just partially differentiate the objective function w.r.t. to the centroids, equate it to zero to find the centroids that minimizes this function? Why do we have to use gradient descent to reach the minimum step by step?


title: How to predict the time gap between two conditions

text: Could anyone suggest the best method to predict the time gap between two events? For example, given that diabetics are at higher risk of developing hypertension, I would like to predict the time gap until patients develop hypertension after being diagnosed with diabetes. I will have both time dependent and independent covariates. Does this require a longitudinal approach?


title: Deep belief nets and PCA

text: I have been using successfully DBNs for classification for MNIST and some other tasks. I was thinking of doing dimensionality reduction before training, in order to boost performance. However, when I tried out using the output of PCA to train a deep belief network the results were quite bad. I understand that DBNs learn features from the image and that PCA might project the data into a space where things do not look as similar anymore, but I was not expecting such a bad result. Any ideas why this could be? I was especially suprised because this paper has succesfully used it: paper here Thank you!


title: Low classification accuracy for statistically different features

text: I am a little bit confused about the machine learning outcomes on my dataset. I would be grateful if anyone can enlighten me on this: When features are statistically different between two groups, shouldn't the machine learning able to predict the right class with a reasonable accuracy? The data consists of 9 features where 8 of them were significantly different between the groups (two tailed independent ttest). Total observations were 71 (30 and 41). However, the machine learning were only able to produce around 65% classification accuracy (10 fold cross validation). I have used SVM with RBF kernel, kNN, Random Forest, and Naive Bayes. My question is, is it normal? Is it possible to get higher classification accuracy when the features were not significantly different between the groups? Thanks.


title: A case of understanding customer behavior

text: Suppose I have a big online company, and many of my customers churned (i.e. they were paying, and then stopped). My goal is to understand why each of them churned. First I identify the complete set of reasons for churning, $H_1,\ldots,H_n$. E.g. "the website is confusing", "the customer lost interest", etc. Then by observation I identify a set of customer behaviors which are correlated with the churn rate, $O_1,\ldots,O_m$. Each of them is directly measurable, e.g. the number of days he logged in, or the amount of time he spent at my website the day he signed up. Assume that computing $\Pr[O_j=\cdot]$ is easy, and that all of my customers are mutually independent. I know that each $H_i$ is related to a subset of the $O_j$'s. E.g. if a churned customer logged in $d$ days, then the probability that he lost interest is $f(d)$, with $f$ a decreasing function. (One of my problems is accurately computing this function $f$.) My goal is to compute, for each churned customer, and for each $i$, the probability that he churned because of $H_i$, i.e. $$ \Pr[H_i] = \sum_j \Pr[H_i\ |\ O_j=x_j] \Pr[O_j=x_j] $$ where $x_j$ is the value of $O_j$ for the given customer. Now, suppose I send a survey to all churned customers, asking for their reasons, but certainly only a few of them answer it. I want then to update all the functions $f_{ij}(\cdot) = \Pr[H_i\ |\ O_j=\cdot]$ based on the new evidence, so that I can compute more accurately $\Pr[H_i]$ for all the churned customers which haven't answered the survey. My question is how can I do this exactly? and how can I get "good" estimates for the original $f_{ij}$'s? Is my approach correct? EDIT: Now I realize that the formula I gave for $\Pr[H_i]$ is obviously wrong. The $O_j$'s don't form a partition. I think I should compute $\Pr[H_i\ |\ (O_1,\ldots,O_m)=(x_1,\ldots,x_m)]$ directly, since the $O_j$'s are not even pairwise independent.


title: Tolerance for pseudoinverse

text: In calculating the pseudoinverse of a matrix $A$, of size (m,n), I need to choose a tolerance threshold for the eigenvalues. I'm trying to understand how I should pick this. Matlab default is to use the $norm(A)*max(m,n)*eps$, where the norm is norm-2 and the eps is the machine epsilon. Any ideas as to why this is the default value? Moreover, there are cases in which $ A $ is only interesting in comparison to other values. For example, I'd like to calculate the eigenvalues of $C + BB^TA^{-0.5}$ where A,B,C are matrices, but B can have a very large second dimension. In this example, how should we choose the tolerance?


title: Dealing with grouped / rounded data

text: I have a dataset that includes variables about customer income levels. The income was collected in binned fashion (`Which range describes your income? 0-25k, 25k-50k,...`). My question is how best use this for modeling using `glmnet` and `gbm` packages in `R`. I have looked at the `grouped` packaged in `R` but it seems to do everything (coursening and regression) for you. Is there a package that converts binned data back to continuous data for use in with other algos? EDIT: The current method I'm using is to convert them to the mid-point of the range (`0-25k -> 12500`), then using an `ifelse()` stmt to code a few variables to convey the fact that there is a relationship between the levels.               incOver25k <- ifelse(df1$income >= 25000,1,0)         incOver50k <- ifelse(df1$income >= 50000,1,0)      Then use these flags instead of using `model.matrix()`. Was curious if there were any better ideas.


title: Fitting a Poisson distribution from missing observations

text: I am interested in fitting a Poisson/negative binomial distribution to estimate the number of times a phenomenon happens within a period, let's just say 10 years. I can count the events from monthly reports, but unfortunately, there are reports missing. So for one sample, I might have 120 observational slots, but for some others I might have 30. The event can happen if it is not observed. The missing slots pattern is random (ie not correlated between samples), and it can result in anything from a nearly complete observational record to a very decimated one. How can I cope with this?


title: R latent class multinomial logit model

text: I am using the `flexmix` package to estimate latent class multinomial logit models in R. In choice theory, there can be variables associated with the alternative (generic) or that vary with the agent (alternative-specific). The `nnet` package that underlies `FLXMRmultinom` can't accommodate generic variables. So far, I haven't seen that the `FLXMRcondlogit` can handle alternative-specific ones. The most flexible package for MNL models is `mlogit`. Has anyone seen an implementation of this for `flexmix`?


title: Decision tree and missing values

text: I need to predict with a simple model if a bank customer is still active. I have some information like the date of the last meeting with the customer, but some of them didn't have any meeting with their advisor.               Customer_Id |  Last_meeting_date (in days from now) | Gender |...     A                115                                   F      ...     B                NA                                    F      ...      Is it possible to use a special value/flag like 1000 to "indicate" to the tree that we don't have any date ? In my mind, the decision tree isn't linear like regression, so it can make an exception, like if Last_meeting_date=1000 then... but I'm not sure of it.


title: Binary logistic regression on SPSS 20 - output seems inconsistent and illogical

text: I have one categorical dichotomous dependent variable (yes/no - retention of newly learned words) and two independent variables - both are categorical. One independent variable is "delay between exposure and test" and has two levels (short time delay/week time delay). The other IV is "word type" and relates to the types of words children were exposed to - object label, colour label etc. I have run a binary logistic regression on my data (440 participants - 44 per condition of word type and time delay) using the "Enter" Method and defined both IVs as categorical using the categorical tab. The initial outputs - Block 0 "Variables not in the Equation", Block 1 Omnibus Tests of Model Coeffcients and the Hosmer and Lemeshow Test all suggest that the model is significant and the Block 1 Classification Table suggests there's been an 11% improvement in predicting the DV (when compared to the Block 0 classification table). Yet the Block 1 "Variables in the Equation" Table shows no main effect of either of the IVs nor an Interaction effect - and not by a long way (time delay p=0.37; word type p=0.37; Word type*Time Delay p=0.401. I don't know how to interpret this as all the examples I've seen in my book and online show at least one significant main effect when the model is significant. It also seems somewhat illogical for SPSS to report that the model is significant but that neither of the IVs are having a significant effect! Also, my data are quite clear and it seems very likely that there is a main effect of time delay. Retention of all word types except one fall from around 65% in the short delay to around 33% after one week. I have tried changing the indicator from first to last (ie from short delay to the week delay) for the time delay variable and concerningly the results change - I thought when the IV was dichotomous it didn't matter which of the categories was used as the control group (indicator). Even more of a concern, it produces a main effect of word type rather than time delay. I have used dummy coding and effectively treated each of my word type categories as different IVs and this produces completely different results. Again - this is a concern as I would expect very similar results since it should effectively be running the same data. Can anyone help?


title: Laplace approximation of the likelihood Bayesian

text: I need help with the following question: Consider $m$ observations $(y_1; n_1); ... ; (y_m; n_m)$, where $y_i \sim Bin(n_i; θ_i)$ are binomial variables. Assume that $θ_i \sim w_1Beta(α_1; β_1) + w_2Beta(α_2; β_2)$ are mixture from two Beta distributions $(w_1 + w_2 = 1)$.   * Derive a Laplace approximation of the likelihood and each mixture component of the prior.    * Derive the empirical Bayes likelihood of the data by integrating out i using the Laplace approximation, and leave the hyper-parameters $(w_j ; α_j ; β_j)$ (j = 1; 2)    * Derive the EM algorithm to estimate the hyperparameters (you may use also mixture of Gaussian prior to approximate mixture of Beta prior)


title: Extremely fast regularized linear model libraries

text: I have a huge dataset (millions of rows, thousands of columns) and `glmnet` Lasso regularised regression is too slow. I am wondering what other libraries there are that try to implement regularised linear model estimation extremely efficiently? I don't necessarily need a distributed solution (it can estimate on one CPU), and it is fine for me to load all data to RAM (so an `R` solution may be fine, if one exists that's faster somehow than `glmnet`). It also does not necessarily need to be Lasso or Elastic Net. I can use java, python or R. My data is pretty messy, some columns are sparse, some have heavy autocorrelation, etc.


title: Possible problem for Network Analysis

text: I have a data set and I want to examine an hypothesis in there and probably Network Analysis should prove or reject my theory. I have a list of products and a group of people who give the product to each other. Anyone can sell it to another in any price he want (hopefully someone will buy it). So, imagine a chain of products transfers between the people. Except the sell price, one will also take a percentage of the sale only of the next one who will re-sell the product. And also, if I sell for the first time a product, I will take a percentage forever for every sale. Example: A -> B -> C Person A will take the money from sale to Person B. Person B will take the money of the sale to person C and Person A will take a percentage of this sale as the previous owner and the creator of the chain. If person C sell the product to someone else, only person B will take the previous percentage, but person A will take the percentage of the creator. I believe that there are people that try to exploit the system and create small or bigger circles to take advantage of these percentages of next sales to "hide" themselves. Is there any way to identify those people?


title: Very large theta values using glm.nb in R - alternative approaches?

text: While analysing the effect of environmental data on the activity of an animal species (the latter given as count data) I am fitting negative binomial GLMs with one predictor using the MASS library in R. Unfortunatley, the data set is very small (n=7 to 9). In some cases, the theta value in glm.nb gets very large (accompanied with the warning "iteration limit reached"), possibly indicating that there's no overdispersion and a poisson GLM might be a better choice. Using a poisson GLM, however, a residual deviance of e.g. 150 on 7 degrees of freedom indicates that there actually is overdispersion - or did I miss something? Using a quasi-poisson GLM works, but I would like to retain ML-based measures such as AIC and Vuong test for model comparison. Any suggestions for alternative approaches are greatly appreciated!


title: Sparse Autoencoder [Hyper]parameters

text: I have just started using the autoencoder package in R. http://cran.r-project.org/web/packages/autoencoder/index.html Inputs to the autoencode() function include lambda, beta, rho and epsilon. What are the bounds for these values? Do they vary for each activation function? Are these parameters called "hyperparameters"? Assuming a sparse autoencoder, is rho=.01 good for the logistic activation function and rho=-.9 good for the hyperbolic tangent activation function? Why does the manual set epsilon to .001? If I remember correctly, "Efficient Backpropagation" by LeCun recommends starting values which are not so close to zero. How much does a "good" value for beta matter? Is there a "rule of thumb" for choosing the number of nuerons in the hidden layer? For example, if the input layers has N nodes, is it reasonable to have 2N nuerons in the in the hidden layer? Can you recommend some literature on the practical use of autoencoders?


title: Diallel Analysis

text: I apologise for the mix-up. Am trying to do a diallel analysis involving three parents in poultry. I tried using SAS05 program but it didnt work, i guess it was designed for diallel analysis involving 4-12 parents. I will appreciate any advice or direction on how to perform diallel analysis with 3 parents using either R or SAS. Thank you


title: do information entropy probabilities have to sum to one?

text: My understanding of information entropy is that it requires the input probabilities to sum to 1. So, for a sequence a,a,b,b you then have -([1/2 log2 1/2] + [1/2 log2 1/2]) = 1 Are there versions of information entropy that don't require probabilities to sum to 1? Or, is there a way to measure entropy that is also sensitive to the quantity of items, not only their probability? Or, is there an accepted way to derive a form of 'non-normalised' information entropy that somehow takes into account that the longer the information stream is, the more likely you'll come across various arrangements of information? e.g. Not that this is accurate, but to convey the question. Lets say you can compute a non-normalised entropy for the same sequence a,a,b,b as such: -([1/2 log2 1/2] + [1/2 log2 1/2] + [1/2 log2 1/2] + [1/2 log2 1/2]) = 2 Alternately, can you sum the information content over a string of information?   * For a,a,b,b you have four items at 1 bit of surprise each, therefore 4 total bits.   * For a,a,a,a,a,a,a,a,a,b you have 10 items at 0.469 bits average surprise, therefore 4.69 total bits?


title: Dynamic Time Warping for irregular time series

text: I have been reading a lot about Dynamic Time Warping (DTW) lately. I am very surprised that there is no literature at all on the application of DTW to irregular time series, or at least I could not find it. Could anybody give me a reference to something related to that issue, or maybe even an implementation of it?


title: Model diagnostics for a glmmPQL in R mixed-effects model

text: Several texts (both online and published books) have been reviewed prior to asking this. What diagnostics are accepted as best practise for a generalised linear mixed- effects model fitted in R using glmmPQL. I am modelling mortality incidence (count data) using repeated measures, longitudinal survival data. The response is either 0 (alive) or 1 (death). Therefore when plotting residuals, a random scatter is difficult to observe as the residuals cluster around the observed values of 0 or 1. How can I plot residuals that are easy to interpret when fitting a glmmPQL model? Is there a way to get one residual per unique subject as opposed to 1 residual per observation of a given subject? What are the other goodness-of-fit measures/tests to apply to a glmmPQL model. E.g. it appears that AIC and BIC can not be obtained from a glmmPQL. I would greatly appreciate any help, even if direction to certain texts. Thanks


title: How do you predict a continuous value from many booleans & a continuous value?

text: Hello: I am a computer science student working as a research assistant in an undergrad IR lab, feeling spectacularly out of my element. Given an input of a single continuous value and a vector of several dozen boolean values, I must estimate the expected value of a single continuous value and output it. I have several thousand training examples where I have _both_ the full input and the actual output. I suspect there is some sort of machine learning algorithm that will produce a function mapping the input to the output, but I'm unsure enough of the terminology involved that I'm not sure how best to seek these resources. What I am looking for reminds me of classic IR classifiers such as the naive Bayesian classifier, but classifiers give you the probability of a sample belonging to a discrete class, not an expected continuous variable. Is this a form of regression? If so, what type? Does anybody have any insight? Any webpages to help fill in the gaps in my knowledge? Or even any helpful search terms? Also, if this is an inappropriate question for stats.stackexchange.com, is there a more appropriate community based QA site you could direct me to? Thanks.