title: Dissimilarity rating experiment: how to analyze the results

text: ### Study: I simulated some surfaces materials at audio and haptic level, and I asked subjects to evaluate on a 9-point Likert scale the degree of coherence between the two stimuli. For example there are stimuli with metal at auditory level and snow at haptic level, or wood both at auditory and haptic level. So it is like a dissimilarity rating experiment. The experiment had only 12 participants, so for each stimulus I have 12 response (no repeated measures involved). ### Question: How should I analyse the following experiment? ### Initial Thoughts: So far the only way I am thinking to use is ANOVA. ### Update Hello again, I used the cmdscale command in R, but I did not get thecperceptual map I would love to see. I as you an help! Maybe the problem is that I misunderstood if my table is well built for the purpose of the analsys. Let´s summarize my experiment, my goals and the problem: I simulated some surfaces materials at audio and haptic level, and I asked subjects to evaluate on a 9-point Likert scale the degree of coherence between the two stimuli. For example there are trials where at auditory level a metal surface was simulated and at haptic level the snow surface, or wood both at auditory and haptic level. The experiment had only 12 participants, so for each trial I have 12 response (no repeated measures involved). In total there were 36 trials, and each trial was evaluated only once by each of the 12 participants. So each subject provided 36 ratings on the 9 point Likert scale, one for each trial. Basically, these are my data:                WD   MT   SW   GR   SN   DL      WD 7.00 6.50 4.91 4.83 5.50 5.00 MT 7.33 6.91 2.08 3.16 4.25 3.25 SW 2.91 1.75 7.91 6.25 6.83 5.41 GR 2.91 2.66 6.25 6.41 7.25 6.75 SN 4.00 4.00 5.58 6.00 7.00 6.58 DL 3.91 3.08 5.16 6.25 6.50 6.83 On the rows the haptic stimuli and on the columns the auditory stimuli. In each cell there is the average score for each trial (e.g. the trial GR-MT is 2.66, that is the average score given by participants to the trial where the material "gravel" was provided at haptic level and the material "metal" was provided ar auditory level). Now I want to analyze the data in the correct ways, and as said MDS is the best analysis instead of ANOVA as I was thinking. My first goal is to print a perceptual map where to place the pairs of audio- haptic stimuli (e.g. WD-WD, MT-DL, etc.) and see how far are the trials from each other. I used cmdscale in R but I did not get the wanted result. Any suggestion? My second goal would be to find some p-values like I normally get with ANOVA. For example I would like to understand if having the coherent couple of stimuli SW-SW (which means "snow" both at audio and haptic level) produces significant differences n the evaluations rather than the couple SW-MT (which means "snow" at audio and "metal" at haptic level) Again I would like to undestand if there is any statistical difference between all the couples of stimuli corresponding to solid surfaces (like the couples MT-MT, MT-WD, WD-WD, MT-MT) and all the couples where a solid surface and a aggregate surface are presented (like the couples MT-SN, or WD-GR, etc.). ...I want to get as many information as possible from that table. I really thanks anyone who can provide any suggestion or useful information.


title: Should percentages have decimal places?

text: When presenting data using a percentage, is it a good thing to have decimal places, say 2 decimal places instead of round off to whole numbers. For example, instead of 23.43%, you round off to 23% I am looking at this from the perspective of whether the 2 decimal places accuracy will make much difference since we are dealing with percentage and not raw data value.


title: Dimension reduction for sparse matrix for clustering

text: I'm looking for a Sparse matrix dimension reduction. I already used some feature selection methods like PCA but it doesn't give me good results. I want to apply mixture models for clustering my data. I appreciate any comment or answer .


title: Visualization of binned frequency distribution in R

text: I am trying to visually represent the distribution of last names in the US. Specifically, I am trying to show that the distribution is such that the most common names (say the top 50) are very common but that it drops off quickly after that. The conclusion I hope to support is that it's not that meaningful to differentiate between common and less common among names that are not in the top X because they are all a very small portion of the population. I have the observed frequency of all last names appearing more than 100 times in the 2000 census.               ## Data from the US Census, extracted and CSV re-hosted      ## http://www.census.gov/genealogy/www/data/2000surnames/names.zip     names <- read.csv("http://samswift.org/files/app_c.csv")      My intuition was to bin the ranked list by groups of 50. Most common names 1-50, 51-100, ...               sum50 <-  tapply(names$count, (seq_along(names$count)-1) %/% 50, sum)      So we now have the sum population of people with a top 50 name, a second 50 name, etc. I was imaging a plot like this ![desired plot](http://i.stack.imgur.com/Nefs2.jpg) where the x-axis is the ordered factor of bins (1-50, 51-100 ..) and the y-axis is the sum population in that bin. I think it's important that bar widths scale with the y variable too so that the area of the square conveys the mass of the population. So, two part question really (although I think that's frowned upon)   1. How might I generate this plot in R with the provided data. I generally use ggplot2, but I am not wed to it. I tried using geom_bar and trying to set the width, but I failed to generate anything even a bit functional.   2. Do you have a better thought on how to visualize the assertion I'm making, or disagree with the assertion entirely?


title: Descriptive statistics for non normal data

text: When I use an histogram to plot my data in Stata almost all variables cluster around one value, _i.e._ these are either at 0 or 1 (a huge bar at one or 0). What does it show and how should I correct it? I am working on panel data.


title: The biggest difference between medical statistics and biological statistics

text: I am turning my work from medical statistics into biology statistics (fishery). I am wondering if there is a big differences in methodology or focus between these two fields? In medical field, we do a lot clinical trials or patient cohort studies, and we try to find associations by applying different types of significance tests and regressions. In the case of biology, there is very limited control of the sampling method. I have a feeling that the statistician's main focus would be to find the best way to sampling the data. Since I am really fresh in this field, can someone give me some ideas about the main focus, common methodologies in biology statistics? Is it really different from medical statistics? Do you have any nice books to recommend? Thanks


title: What does this formula (to derive the annualized volatility for the VaR) mean?

text: I'm faced with the formula shown in the image below, which I just don't understand, in part because I've no grounding in stats, and in part because I don't even understand the notation: ![enter image description here](http://i.stack.imgur.com/qrITP.png) What's going on here? Is this showing in the first line how to compute the VaR, and then in the second line how to derive the annualized volatility from some of the variables used in the first computation? I don't even understand the "T time intervals of 1/m years" part. What's an "mth" of a year? And what does "(sigma) 1/m" mean? Is that the volatility for _one_ time period? If so, then don't we need to aggregate the values for all the discrete time periods somehow? I'd expect to see a "sum from 1 to m" somewhere... I'm totally confused, and any help would be very gratefully received.


title: When to use non-parametric regression?

text: I am using PROC GLM in SAS to fit a regression equation of the following form $$ Y = b_0 + b_1X_1 + b_2X_2 + b_3X_3 + b_4t $$ The QQ plot of the resulting redsiduals indicate deviation from normality. Any transformation of $Y$ is not useful in making the residuals normal. At this point, can I safely switch to non-parametric methods such as PROC LOESS. I have already used PROC LOESS, and the fit looks better than PROC GLM. But I do not have much knowledge in non-parametric regression. I do not know when to choose non-parametric regression over parametric regression. Can someone help me with this? I will go ahead and add another question. Following are the description of my variables in the model. $$ Y =\text{cost of medical care}\\\ X_1 =\text{number of injections}\\\ X_2 =\text{number of surgeries}\\\ X_3 =\text{number of physical therapies}\\\ t =\text{time} $$ Sometimes I get negative predicted cost. This does not make sense. How can I address this issue?


title: Harmonic mean with zero value

text: How does harmonic mean handle zero values? what would the harmonic mean of {3, 4, 5, 0} be since $1/0=\infty$?


title: To what extent can we call a Geometric Distribution a Geometric Density

text: In some papers, for example in "The Geometric Density with Unknown Location Parameter" by Klotz, a Geometric Distribution is called a Geometric Density. For me, this claim looks erroneous, however Klotz is a serious statistician and a professor in the field. My question is, to what extend is it legitimate to call a Geometric Distribution a Geometric Density?


title: Pattern detection in a time series

text: Lets say I am trying to find patterns in pricing of a commodity such as corn and I have monthly time series data on the demand, the temperature, cost of water, and the current market supply. What methods would I use to detect patterns such as: *Corn will cost less during harvesting in October and November when temperatures are around X. *Corn prices are less when the supply goes up on the month of the Corn festival. *Last year corn prices are high because there is a drought and the cost of water and the temperature is high in the summer.  What would I use for these types of patterns in a time series? Usually you would use a correlation matrix for linear data. But what do you use for seasonal and stochastic data where there are patterns like these exist but you do not know what they are? Are these types of cases for cointegration and vector autoregression (my guess)?


title: E-M, is there an intuitive explanation?

text: The E-M procedure appears, to the uninitiated, as more or less black magic. Estimate parameters of an HMM (for example) using supervised data. Then decode untagged data, using forward-backward to 'count' events as if the data were tagged, more or less. Why does this make the model better? I do know something about the math, but I keep wishing for some sort of mental picture of it.


title: Error propagation on correlated Poisson samples

text: If I have a set of $n$ tuples ($b_1,b_2,b_3,b_4,x,y,z$), where $b_1, b_2, b_3$ and $b_4$ are counters ruled by $P_1, P_2, P_3$ and $P_4$ Poisson distributions (error on the counters are $\sqrt{b_1}$, etc.) and $x, y$ and $z$ are ruled by real valued distributions. From it I select a subset of tuples $m$ according to a condition on $x,y$ and $z$, e.g. $x > P_I$. That subset will have certain $b_1', b_2', b_3'$ and $b_4'$ values for the counters. If e.g. I compute the $b_1'\div b_1$, what's the error on that ratio?


title: Something like E-M for discriminative models?

text: E-M provides a way to improve the estimation of a generative model with unannotated data. Is there anything out there that works the same way for discriminative models (e.g. perceptrons)? For example, consider averaged perceptron tagger. It would be handy to be able to throw the entire Gigaword through some process of unsupervised model improvement. **EDIT** : So, I was pleasantly surprised to note that this site has the ambition of dealing with machine learning, but I'm learning by experiment what vocabulary is generic and what is very domain-specific. Apologies. Consider a sequence classification problem, like part-of-speech tagging or named entity extraction. You can train a generative model (e.g. an HMM). Thats's a probability model, and you can apply E-M. However, the number of states grows prohibitive if you want to look at many features, and so the fashion tends toward things like CRFs (batch) or Perceptron (online). For example, this paper talks about unsupervised learning in for a perceptron POS tagger, but the details are that they add the output of several pre- existing taggers to the training set of their model.


title: Conflicting Normality Tests

text: I have been trying to establish whether I could model my data, a set of 90 observations, using the normal distribution. I have tried the Shapiro-Wilk and the Anderson-Darling tests but each came back with a different result. **Anderson-Darling normality test** data: x A = 0.6994, p-value = 0.06555 **Shapiro-Wilk normality test** data: x W = 0.9645, p-value = 0.0154 If we are to use the strong 1% significance level, the null of normality cannot be rejected in either case. Using the 5% size, however, we see that the SW test rejects the null while the AD one narrowly accepts it. The results are conflicting and I do not know on which to place my trust. If I knew the power of each test, I could make up my mind so is there perhaps a way to obtain the power functions in R? The Normal distribution is easy to work with and it would make sense in my case seeing that my observations are test scores but I do not want to make farfetched assumptions. I also have to say that personally, I am not convinced that the underlying distribution is normal. Please take a look at the histogram below, upon which I have fitted a normal distribution with the same mean and standard deviation. ![enter image description here](http://i.stack.imgur.com/pkncK.jpg) It does not seem to be a good match. Below is the normal qq plot: ![enter image description here](http://i.stack.imgur.com/7h5Tc.jpg) There is a hint of linearity but there are also a lot of outliers. As far as I know, the Shapiro-Wilk test is based on order statistics and it does not support the normality hypothesis. If we were to reject the normality of my data, then what would be an acceptable distribution here? Note that we are looking for a negatively skewed one. All suggestions are welcome. Thank you. **EDIT** : This is my vector of observations in ascending order, i.e. the 89 order statistics, in case anyone wants to test it further: > sort(x) 19 32 37 37 37 42 42 45 45 45 46 46 46 47 47 48 48 50 51 54 55 55 55 55 55 55 55 56 56 56 56 57 57 58 58 58 60 60 60 60 61 61 62 62 62 62 62 62 62 62 62 62 63 63 63 65 65 65 66 66 67 67 67 67 68 68 68 68 70 70 71 71 71 72 72 75 75 76 76 77 77 77 77 77 78 78 78 80 81


title: Assessing questionnaire reliability using test-retest measures

text: I am an undergraduate student with not much knowledge regarding statistics. I was hoping you would be able to help me by answering a question. I am currently checking the reliability and validity of the questionnaire. For the reliability of the questionnaire, I was going to carry out a "test retest" with two students giving them the same questionnaire with a 2 week lag to see if they answer in the same manner. I would like to know if Cohen’s kappa coefficient can be used for my questionnaire. I have done some research and found the equation to be: Kappa = Proportion of observed agreement — chance agreement /1 - chance agreement In my case, would the 'chance agreement' be still 0.5 seeing that, in my self- made questionnaire there are many options to each question and it is not a simple yes or no question. If this is not what I should be using, can you suggest any other method to carry out the reliability of the questionnaire?


title: How to interpret and compare models in Cox regression?

text: I am trying to interpret the results of a Cox regression; I am doing a PhD in medicine. I love statistics but my question is still pretty basic, I think, and I did not find an answer in previous threads. I have to compare different models (just a couple of predictors in each; predictors are different, but sometimes the same predictor appears in different models; say, A+B, A+C, A+D, B+C) toward the same time-to-event variable. How do I choose the “best” model? I am studying the underlying statistical principles, but I still don’t get whether I have to look at   * Which model has the highest log likelihood;   * Which model has the best p-value of the likelihood ratio test (LR chi2);   * Which model has all p-values of the HR coefficients (beta’s) of the covariates significant;   * Or any combination of the above (for instance, only consider the models where both the LR and the beta-coefficients of all covariates are significant, and among them choose the one with the highest log likelihood). Or is there another statistical technique to do it? I am using Stata 11.0. I understand that the significance of the LR and the significance of the beta- coefficients test different things, but still I need to select the model with, say, the “best predictive ability” or the “strongest association”. Thank you in advance for your help! Luca


title: Predicting time series with NNs: should the data set be shuffled?

text: Suppose I'm trying to predict time series with a neural network. The data set is created from a single column of temporal data, where the inputs of each pattern are `[t-n, t-n+1, ... , t]`, `t` being the time step and `n` the embedding size, and `[t+1]` being the target (predicting the "next step" of the series). Here is the question: if I use such a data set for NN training, should I shuffle it? Shuffling is usually done to avoid overfitting by preventing the NN from memorizing the pattern order. However, in case of time series, could pattern order be a form of information that I shouldn't deprive the NN of? The issue becomes even more severe if recurrent NNs are concerned, where every pattern recurrently receives information about the previous pattern: would a RNN benefit if patterns were in order?


title: Plot regression adjusted for covariate

text: I would like to plot values from a linear regression adjusted for a covariate. Please consider the following dummy data in which y is predicted by x and the covariate a. There is an interaction between x and a.               # generate some dummy data     x1 <- rnorm(20,10,5)     x2 <- rnorm(20,20,5)          a1 <- rnorm(20,10,5)     a2 <- rnorm(20,20,5)          y1 <- x1*a1+rnorm(20,5,3)     y2 <- x2*a2+rnorm(20,5,3)          x <- c(x1,x2)     y <- c(y1,y2)     a <- c(a1,a2)          # create model with & without covariate     m1 <- lm(y~x*a)     m2 <- lm(y~a+a:x)     m3 <- lm(m2$residuals~x)          summary(m1)     summary(m2)     summary(m3)          # plot data     par(mfrow=c(1,2))     plot(x,y)     abline(lm(y~x, subset=c(1:20)))     abline(lm(y~x, subset=c(21:40)))     plot(x, m2$residuals)     abline(m3)        * Is this an appropriate way to visualize the data?    * Any other ideas how to visualize a main effect when interactions are present?   * Why is in model m1 a main effect of x but when I analyze the residuals of model m2 in model m3 there is no effect of x?    * Is that as inclusion of different predictors will always change the significance of other predictors?


title: Network/structure learning

text: Given a data set $\mathbf{X}\in\mathbb{R}^{n\times p}$, where $n$ is the number of samples (observations) and $p$ is the number of features, I would like to know what kind of methods exist for learning the hidden network of relationship among the features and detecting from this network the relevant groups. I am a bit confused about the output of probabilistic graphical models and structured-sparsity regularization methods. Is it possible to use probabilistic graphical models (Bayesian methods or Graphical Lasso) to learn a network from a data set? Is it possible to use structured-sparsity regularization methods (group Lasso, graph Lasso, etc.) to uncover relevant structures (communities) from the learned network? Are these two approaches complementary? Could you suggest me any relevant papers?


title: Building a customer cancellation predictor

text: I'm trying to build a cancellation predictor for telecom data. I am using both static (i.e. location, device, number of complaints, etc.) and temporal (i.e. time-series usage) data. The response variable is whether or not they cancelled within the first 2 months of activation. My question is related to the response variable. In addition to whether they have cancelled, I also have access to a short text field that the customer has filled out where they explain why they have cancelled. In some cases, the cancellation reason seems legitimately unpredictable (eg. "My phone got stolen" or "I am merging with my wife's account" or "moving out of the country" or "I was just testing the service"). I was wondering if I should remove these instances from the training set. It seems like no matter how good the prediction model is, it can never predict a phone getting stolen for example. So it seems like if I leave these in, it would only hurt the model's performance. Is this a correct assumption?


title: Whether to exclude a case from two analyses when the case is only an outlier in one?

text: I have a question regarding the elimination of participant data for two statistical analyses on the same sample. Here is a bit of background information: I collected data from one sample (n = 181 college students) and administered two questionnaires to them (one on spirituality and one on negative emotional states). I performed an exploratory factor analysis on the data from the spirituality questionnaire, and I want to compare the factors derived from the EFA to the negative emotional states using a hierarchical multiple linear regression. My main question is the following: When I analyzed the data with a hierarchical multiple linear regression, SPSS identified (via case diagnostics) participants whose data could be eliminated due to outlier scores. If I exclude these cases from the analysis, then do I have to remove them from the results of the EFA as well? (Both of these analyses would be featured in the results section of the same experiment.)


title: Values p,d;q in ARIMA

text: If the process is: $(1-0.3L)(1-L)X_t=(1+0.5L)^2\epsilon_t$ what are the values of p,d and q in ARIMA (p,d,q)? Thank you very much in advance.


title: computing slope on heavy tailed lognormal distribution

text: I have a heavy tailed logNormal distribution and i want to know if it's possible to make a linear regression in order to compute R² and slope for this data. I make my linear regression on computed values log10(rank) and log10(population) (cf it's a rank size plot) ![enter image description here](http://i.stack.imgur.com/TQ01x.png) I don't know if it's a good method, especially in the case of heavy tailed distribution. Do you know better/robust method to compute slope of lognormal distribution?


title: Three versions of discriminant analysis: differences and how to use them

text: Can anybody explain differences and give specific examples how to use these three analyses?   * LDA - Linear Discriminant Analysis   * FDA - Fisher's Discriminant Analysis   * QDA - Quadratic Discriminant Analysis I searched everywhere, but couldn't find real examples with real values to see how these analyses are used and data calculated, only lots of formulas which are hard to understand without any real examples. As I tried to understand it was hard to distinguish which equations/formulas belonged to LDA and which to FDA. For example let's say there is such data:               x1 x2 class     1  2  a     1  3  a     2  3  a     3  3  a     1  0  b     2  1  b     2  2  b      And let's say some testing data:               x1 x2     2  4     3  5     3  6      So how to use such data with all these three approaches? It would be best to see how to calculate everything by hand, not using some math package which calculates everything behind the scenes. P.S. I only found this tutorial: http://people.revoledu.com/kardi/tutorial/LDA/LDA.html#LDA. It shows how to use LDA.


title: What are the rules for picking a distribution for modelling a data set?

text: I have a data set of 100,000+ event times. For this study an event begins (at time 0) and can run for an indeterminate amount of time. The bulk of events require a couple of hours to complete but there are also events that can take tens of hours to complete. I created a histogram of the times and it resembles a positively skewed distribution. My question is: how can I determine a distribution type to model this data? For example, the data looks log-normal, but a histogram of the log of the data is not a normal distribution (it is a negatively skewed distribution). I spent some time looking at the Gamma distribution but besides visually resembling my data I don't know if Gamma is an appropriate selection. I also had a look at the type-2 Gumbel distribution but only because it is asymmetric and semi- infinite, like my data. What are the rules for picking a distribution to model a data set?


title: Should the likelihood function be increasing in every step of the EM algorithm?

text: Should the likelihood be increasing in every step of the EM algorithm? I wrote an EM algorithm recently and the number it arrived does not seem to be the maximum. I know this because I used the optim function on R on the same problem and it arrived parameters which gave a bigger log-likelihood value. But even when I started at this maximum, I run my EM algorithm, it seem to converge to the other value I initially derived from my EM algorithm. What is the most likely explanation for this?


title: In logistic regression, does the lack of significance of the parameter estimates in a test sample indicate overfitting?

text: I am trying to build a logistic regression model where I have a dependent variable $y$ and independent variables $x_1$, $x_2$... $x_n$. $y$ can take only two values - 0 or 1. My original modelling dataset has 100,000 observations - which I have divided into two samples - a training sample with 80,000 observations and a test sample of 20,000 observations. The samples were created randomly, maintaining the same proportion of 1 to 0 in both the samples (using the SURVEYSELECT procedure in SAS and $y$ as STRATA). Let us assume that the percentage of observations with the value 1 for $y$ in both the samples is 10%. I use the LOGISTIC procedure in SAS on the training sample to arrive at variables $x_1$ - $x_5$ which constitute my final model. The p-values associated with the Wald Chi Square are all <0.0001, which indicate that they are significant at the 99.99% confidence level. However, when I run the LOGISTIC procedure on the test sample, using only $x_1$ - $x_5$ as independent variables, one of them say $x_4$ is no longer significant in the regression outputs - that is, the associated p-value is 0.6. Does this mean my sampling is not proper? Or the model which I obtain from the training sample 'overfits' the data? Or both? Should I be worried about this and ensure that the final variables I choose are significant in both the samples or this is not an issue in general?


title: Relaxing the parallel lines assumption in a proportional odds model

text: I tried to specify a partial proportional odds regression in STATA using the `gologit2` command. However, `gologit2` runs extraordinarily slow in my dataset (108k observations of 9 vars). For example,               ologit PfSt age if bund == 1      runs in less than 1 sec, while               gologit2 PfSt age if bund == 1, pl      which is essentially the same, as the `pl` option enforces parallel lines in all vars, i.e., the `ologit` specification, takes about a minute on my computer. Calculating the fully specified model is possible with `ologit` and `mlogit` (~ 1 min), but not with `gologit2`, even if I specify for which variables parallel lines should be assumed. The regression command is               ologit PfSt age age_sq gender i.bundesland tt tt2 [fw=size], cluster(bundesland)      where `PfSt` ranges from 0 to 7 and `tt` capture time trends (`tt = year − 1996` and `tt2 = -1/tt +1`, a “phasing out” trend) Now, considering that relaxing the parallel lines assumption for some variables would be identical to adding interaction terms to my model specification, I could change my regression command to `ologit PfSt age* gender i.bundesland#i.PfSt tt*, cluster(bundesland)` would yield me different coefficients per level of `PfSt` for the one variable for which I would relax the parallel lines assumption. **Is that correct?** _However_ , the estimation process **fails to converge** , as now (1200th iteration) negative log pseudoelikelihood ratios rise (i.e., drop in absolute values, now <<0.0001), but they are `(non concave)`. I wonder **why is this the case?** , as both `ologit` and `mlogit` did produce results after a few iterations.


title: Fitting data/determining parameter values for custom (biPareto) distribution

text: I have a dataset containing information about a bunch of wireless devices. Specifically, the number of other wireless devices each device encounters. Based on other researcher's prior work in the same area, I have reason to believe my distribution may potentially be well fitted by the not so well known biPareto distribution, whose CCDF is defined as: $1 - F(x) = (x/k)^{-\alpha}(\frac{x+kb}{k + kb})^{\alpha-\beta}$, for $ x > k$ and simply $1 - F(x) = 1$, for $x \le k$. Assume I already know $k = 1$. This leaves a function to fit with the three parameters $\alpha$,$b$,$\beta$, with the restriction that each parameter is $> 0$. What I am unsure about is how one would go about determining the most suitable fitting parameter values, ideally in `R`. I am quite the novice when it comes to fitting data, so aside from the syntactical considerations in `R`, I am also interested in any advice regarding what specific methods of fitting would be most appropriate for this particular function. For example:   1. Does this function qualify for linear regression?   2. Is least squares appropriate/applicable here? What about minimum mean squared error (MMSE)?   3. How sensitive is any given fitting method to any "guiding" initial parameter values I give it (iterative methods?), or is the result determined analytically (always produces the same answer in the end?)  If it's at all relevant, I'm ultimately looking to use the Kolmogorov-Smirnov test on the data to see how well it matches the parameterized biPareto distribution.


title: Using an If condition inside SPSS syntax

text: I have a data table where I need to calculate the number of days each user performs an action. The structure is like this.               UserID       Date      TimesActionDone       A        20110601        3       A        20110605        0       A        20110609        10       B        20110701        0       B        20110710        0       B        20110801        1       C        20110801        0      So, the result will be                 UserID     Days         A          2         B          1         C          0      This is because TimesActionDone was 0 for one day for A and two days for B and C is zero. These days should not be counted. How can this be done in SPSS syntax?               AGGREGATE       /OUTFILE=*        /BREAK= UserID Date            /N_BREAK=N      /TimesActionDone =sum(number_of_actions).      This generates the first table. How do I get to the second? I have gone as far as this but doesn't work.               AGGREGATE       /OUTFILE=*        /BREAK= UserID        IF(TimesActionDone > 0  ).      /N_user_days=N.       I have not much knowledge of SPSS syntax and couldn't find much to go with my searching online. I know I can do this by removing all the cases where TimeActionDone is 0 first and then do the calculation but then I would loose information about user C. I also need to have this as 0.


I need a calculator that uses a contigency table of 11x6. Anyone knows where I can find a calculator for my case? Thanks in advance.


title: Hessian of profile likelihood used for standard error estimation

text: This question is motivated by this one. I looked up two sources and this is what I found. A. van der Vaart, Assymptotic Statistics: > It is rarely possible to compute a profile likelihood explicitly, but its > numerical evaluation is often feasible. Then the profile likelihood may > serve to reduce the dimension of the likelihood function. Profile likelihood > functions are often used in the same way as (ordinary) likelihood functions > of parametric models. Apart from taking their points of maximum as > estimators $\hat\theta$, the second derivative at $\hat\theta$ is used as an > estimate of minus the inverse of the asymptotic covariance matrix of e. > Recent research appears to validate this practice. J. Wooldridge, Econometric Analysis of Cross Section and Panel Data (the same in both editions): > As a device for studying asymptotic properties, the concentrated objective > function is of limited value because $g(W,\beta)$ generally depends on all > of $W$, in which case the objective function cannot be written as the sum of > independent, identically distributed summands. One setting where equation > (12.89) is a sum of i.i.d. functions occurs when we concentrate out > individual-specific effects from certain nonlinear panel data models. In > addition, the concentrated objective function can be useful for establishing > the equivalence of seemingly different estimation approaches. Wooldridge discusses the problem in wider context of M-estimators, so it applies to maximum likelihood estimators also. So we get two different answers for the same question. The devil in my opinion is in the details. For some models we can use hessian of profile likelihood safely for some models not. Are there any general results which give conditions when can we do that (or cannot)?


title: AICc results in R

text: I used the model:               Model1 <- glm(disease ~ Distance_ue + Pop_density + Alt + Ave.temp + Ave.precip +                              Urban.Rural + (1| Sites), family=gaussian)      I then used the `dredge` function to obtain an AICc results table. Model selection results for models predicting the effect of various predictor variables on `disease` (presence or absence) in birds in 10 research sites. `Site` was used as a random term. Predictor variables were classification of `Urban / Rural` (categorical - either labelled as urban or rural), `distance to the urban edge` (continuous), `population density` (continuous), `altitude` (continuous), `average temperature` (continuous) and `average precipitation` (continuous). The issue I have is that in the results table there is no parameter estimate given for classification of urban or rural which is the categorical variable - it only gives a $\pm$ sign (please see image below for the results table). Is this normal? Why is it doing this? ![AICc results table](http://i.stack.imgur.com/6c48U.png)


title: Predicting binary vectors as a single response

text: My responses consists of vectors of 12 bits - each bit means a particular event in 1 month. I need to predict the probability of an event in all months -> so the predictions should also consist of 12 probabilities. Examples of responses:               0 0 0 0 0 0 0 0 0 0 0 0     0 1 0 0 0 0 0 0 0 0 0 0     0 0 0 0 0 1 0 0 0 0 0 0      Is there a way to model those 12 bits as 1 response? I was thinking about modeling the month of event so the response would be 13 (meaning no event in 12 month),2,6 etc. But then how do I get the probabilities for each month? The reason I want to do this is not to create 12 independent models but 1 model that understands the distances between months ie probability of an event in January should be simillar to February but much different than December. Maybe there is a standard way to model this. I was reading about survival analysis but I don't know if it is the best option. I would be grateful for any hint.


title: How can I make my data fit normal distribution?

text: I have selected a number of words and asked people to rate them on a Likert scale. Words were selected according to two categories; syllable length (mono-, di- and trisyllable and final phoneme (vowel or consonant). These variables interact and are three words for each cross-category (eg. 3 monosyllabic vowel-final words, 3 disyllabic vowel final words, 3 trisyllabic vowel-final words etc). I am now analysing my data and however I cut it I cannot get it to fit normal distribution. I have organised it by broad categories like syllable length or final phoneme, and finer categories like monosyllabic vowel-final words. I have applied logarithmic transformations and all the p-values returned on the Shapiro and Wilcox tests are absolutely tiny. Am I doing something wrong? Or is this just the way the data is?


title: Sequential linear estimation: Is MSE non-increasing?

text: Suppose we estimate $x$ from $y_1=ax+z_1$ and get $x_1$. Then suppose we receive $y_2=bx+z_2$ and reestimate $x$ using both $y_1$ and $y_2$ and get $x_2$. Is the means squared error (MSE) of $x_2$ always lower than MSE of $x_1$. In other words is $$E(x-x_1)^2\geq E(x-x_2)^2?$$ Here $a,b$ are know constants and $z_1$ and $z_2$ are Gaussian errors and the estimation is linear minimum mean squared error. The equations can also be vectors in that case error is the trance of the error covariance matrix. By this wiki it seems to be decreasing. Edit: explanation At time $t=0$ a real value $x$ is generated by a random process $X$. Then at time $t=1$, I receive a noise corrupted and scaled version $y_1$ of $x$ denoted by $y_1=ax+z_1$. From $y_1$ I do linear minimum mean squared error (LMMSE) estimation and estimate $x$ and get $x_1$ and the corresponding LMMSE is $E(x-x_1)^2$. Here $a$ is a constant and is known, $z_1$ is Gaussian noise with known mean and variance. LMMSE estimation is straightforward. Then at time $t=2$, I receive another noise corrupted and scaled version $y_2$ of that same $x$ which is denoted by $y_2=bx+z_2$. Again $b$ is a known constant and $z_2$ is Gaussian noise with known mean and variance. Now I simultaneously use both $y_1$ and $y_2$ to estimate $x$ by LMMSE estimation and I get the value $x_2$ and the corresponding error $E(x-x_2)^2$. My question: is $E(x-x_1)^2\geq E(x-x_2)^2?$. i.e., the time actually does not matter the question is does LMMSE estimation using $n+1$ samples always give a less error than estimating only using $n$ samples out of that $n+1$ samples? Thanks.


title: Problem in finding fractal dimension (correlation dimension)

text: I am trying to implement Takens correlation dimension formula given in Eq (8) in Paper1: Estimation of the dimension of noisy attractor (Estimation of the dimension of a noisy attractor. Schouten JC, Takens F, van den Bleek CM Physical Review E, vol 50, number 3, 1994). The actual derivation and explanation is given in Takens maximum likelihood estimation (F. Takens, On the numerical determination of the dimension of an attractor, in: D. Rand, L.S. Yong (Eds.), Dynamical Systems and Turbulence, Warwick, 1980, Lecture Notes in Mathematics, Vol. 898, Springer, Berlin, 1981.) and summarized below according to my understanding. My time series are noisy with additive white Gaussian noise, `Z = y+noise` Step1 = Perform phase space embedding of the 1D time series. This becomes U. Step2: Find euclidean distances or maximum norm distances of the higher dimensional space. Find a scaling distance, $l_0$ (or $r_0$) called as the maximum scaling distance. Step3: Then find those distances, $l_z$, which are $\le l_0$. This step converts distances into $l_z$ that now belongs to range $[0,1]$. Sometimes, the maximum scaling length $l_0$ is taken as the $l_0$ = average absolute deviation of the time series by the formula $l_0 = 1/{\rm total\ elements}[ {\sum (u_i - {\tt mean}(u)}]$. Step 4: In Takens paper, it is mentioned under Section2, that the distances can be made independent of $l_0$ by dividing all distances $l_z$ by $l_0$. Step5: Take the logarithm of each element so found by Step 4. ${\rm Takens\ dimension} = -(\frac{2}{n(n-1)} \sum [\log (||X_i- X_j||/r_o)]^{-1}$ (Formula for Takens) $n$ = number of samples. But, I followed Eq (8) of paper1. **Q1:** It is not clear to me what should be the $l_0$ - (A) When I took $l_0$ as the average absolute deviation and divided the distances by it, the result of correlation dimension is absolutely incorrect. So, if I follow what is mentioned in the paper, then the value of the dimension is coming out to be incorrect. What should be $l_0$ ? (B) When I normalize the time series by standardization and then I choose $l_0$ as the standard deviation of the distances (since $l_0$ is labelled as maximum scaling "distance"), then the result comes somewhat close but not exactly precise. (C) When I normalized the time series by `normalized_timeseries=(data-min(data))./abs(max(data)-min(data));` and then $l_0$ is selected to be same as in (B), I get correct answer! But, this violates/contradicts what is mentioned in the paper. I really do not understand what should be $l_0$, the upper bound of scaling distances. does not should be the maximum distance, or the maximum standard deviation of the distance vector ? The value of correlation dimension from Takens method for Lorenz should be approx 2.14. I am getting incorrect value. Please help in implementing the above formula, I am not sure if I understood correctly. I have partly implemented the code, upto the calculation of the euclidean distances.               clc     clear all          iter= 1000;     dt=0.02;     [x y z]=lorenzo(iter,dt);  %generate time series from Lorenz     d=3;     l2=iter-d;     Z=y;          [N col] = size(Z);    % Current size of Z     zmean = mean(Z); % Mean of each component          Z=(Z-min(Z))./abs(max(Z)-min(Z));  % normalization by standardisation     l0 = max(std(normalized_Z));                    for k =1 : l2       U(:,k) = normalized_Z(k+d-1:-1:k)'; %phase space embedding     end          for i=1:l2       for j=1:l2         D(i,j) = sqrt(sum((U(:,i) - U(:,j)) .^ 2));   %Euclidean distance       end      end      Function for Lorenz               function [x y z]=lorenzo(iter,dt);     x(1)=2.9; y(1)=-1.3; z(1)=25.9;     for i=2:iter         [x(i) y(i) z(i)]=nextit(x(i-1),y(i-1),z(i-1),dt);     end          function [XO YO ZO]=nextit(X,Y,Z,dt)         r=28; s=8/3; t=10;         x1=X+t*(Y-X)*dt/2;         y1=Y+(X*(r-Z)-Y)*dt/2;         z1=Z+(X*Y-s*Z)*dt/2;         XO=X+t*(y1-x1)*dt;         YO=Y+(x1*(r-z1)-y1)*dt;         ZO=Z+(x1*y1-s*z1)*dt;


title: How do I compare correlation coefficients of the same variables across different groups?

text: I am testing the correlation between Variable A and Variable B by gender and marital status. So I have the coefficients for males and females and single, married and other. I want to compare the correlation coefficient of male with female. Similarly, I want to compare the correlation coefficient of single with marrried and with other. Finally, I want to compare the correlation coefficient of single male with single female, married male with married female and other male and other female. I have these coefficients too. **What is the easiest way of making this comparison?**


title: Differentiate a ARMA(p,q)-GARCH(p,q) or not?

text: I'm confused about the stationarity condition. I'm fitting a ARMA to a time series that might not have a constant mean. When I fit a ARMA model to it the residuals looks stationary. Do I have to differentiate the time series before I fit the ARMA model or not? I'm confused, since when I look at examples in SAS and r, no one seems to diffirentiate the time series when fitting a ARMA- model, even doe the time series exhibits a trend. My time series does not have a clear trend! Thanks!


title: Quantile regression vs. Li's regression: which should I use, and when?

text: Is there a general rule of thumb about when robust regression or quantile regression is preferred in the presence of outliers? For example, I have a dataset where the DV exhibits extreme positive skewness. However, the large cases are actually some of the most interesting observations. When I run OLS, I find a positive relation between the DV and the IV of interest. When I estimate quantile regressions, I find that the positive relation between the DV and IV of interest is strongest in the 85th, 90th, and 95th percentiles (which is where one might expect). It is insignificant and sometimes negative for the rest of the percentiles. However, when I run `rreg` in Stata, it basically gives no weight to the large positive outliers, leading to no relation between the DV and the IV of interest. Which approach (OLS, Quantile, `rreg`) should be reported? Which is most appropriate?


title: Similarity between two sets of random values

text: I have two sets of random floats between [0,1], each one with a different number of elements: $$A = \\{0.3637852, 0.2330702, 0.1683102, 0.2127219, 0.0152532, ..., N_A\\}$$ $$B = \\{0.4541056, 0.7521812, 0.0266602, 0.5099002, 0.3468181, ..., N_B\\}$$ where $N_B > N_A$. I need to asses the similarity between these two sets _disregarding the difference in number of elements_ (ie: the fact that $N_B > N_A$ should not play a role in the similarity assesment) since I'm only interested in **how** the values are scattered between [0,1], not **how many** of them there are in each set. So far I've applied the 1D `kde.test` available in Duong's R package 'ks' (see page 27) which returns a `p-value`, and I've also applied the `python` function scipy.stats.ks_2samp which is a 1D Kolmogorov-Smirnov statistic which returns a "KS statistic" and a "two-tailed p-value". My questions are: 1- Is one of these statistics (KDE's p-value, KS statistic or the two-tailed p-value) recommended for my needs? If so, why? 2- What is the difference between the "KS statistic" and a "two-tailed p-value"? 3- Will the difference in the number of elements in each set affect the outcome of these statistics? If so, how can I avoid that? Is it even possible?


title: Cross-matrix inference and correlations

text: I have two matrices with the same individuals in each.   1. How should I proceed in R to identify the cross-matrix relations?   2. What technique to use to infer the correlated values in the second matrix from the first one? Thank you.


title: Interrupted time-series analysis for panel data

text: I am familiar with using regression with ARIMA errors to model interrupted time-series, in order to estimate the change in magnitude caused by a policy intervention. These models seem to be designed for a single time series, and thus if multiple time-series are analysed a model must be fit separately for each time-series. I am interested in analysing the national impact of a policy intervention, the implementation of which was staggered in time across all (eight Australian) states. I can see three possible analysis approaches here:   1. Fit a separate ARIMA model for each state   2. Attempt to fit an aggregate national model, perhaps with one dummy variable indicating partial implementation and another indicating complete implementation   3. Find a different model that works explictly on panel data. This would hopefully bring some kind of compromise between the no-pooling approach 1 and the complete pooling approach 2. What sort of approach would you recommend here?


title: Best way to interactively guide a remote statistical session, for a computer without admin access?

text: I collaborate with colleagues in other cities using other operating systems. My colleagues are in the public service and so have limited flexibility to configure their (windows) machines. I want to be able to show them what I'm doing and ideally see what they are doing. How have people tackled this problem? We presently use R but other options are possible. Thanks! Andrew


title: Visualization of a nondirected graph with nodes in specific positions

text: I have a variety of brain regions and information about how they're connected. (Not FMRI in this case, but the diffusion information about the white matter tracts connecting gray matter regions.) I'd like to graph this information sensibly to deliver to the scientist who asked me for this information along with a spreadsheet describing each connection. I'd like to represent both the positions of the gray matter regions (I can project them onto the mid-sagittal plane to get 2D positions or do something similar) and the characteristics of the connections. (Probably just 2 or maybe 3, so I imagine color and thickness would cut it.) I'm not very familiar with graph visualization tools -- I started with the graphviz dot format -- but now I'm reading about a multitude of tools. What sort of graph visualization tools ought I to use?


title: Simple distance measure for financial time series

text: I have a large quantity of financial trading systems that I believe are highly duplicative, meaning that I believe a large number of the trading systems are essentially the same thing. I am looking for a way to measure the differences between each of the trading systems. One simple way to do this is correlation of daily returns, although that limits me to linear relationships and a correlation matrix of 100k x 100k, which is not really useful. **What methods can I use to measure differences between the time series and then to cluster them together,** with the eventual goal of removing the trading systems that are highly similar? I'm new to this, so if I'm leaving out relevant information please let me know and I will modify this question. Thanks. Sample data is below:               StratID SystemID Date   Daily Simple Return     1   1   1/25/2011   0.04     1   1   1/26/2011   0.49     1   1   1/27/2011   -0.02     1   1   1/28/2011   0.76     1   1   1/31/2011   0.61     1   2   1/3/2011    1.37     1   2   1/4/2011    -0.02     1   2   1/5/2011    -0.52     1   2   1/6/2011    0.16     1   2   1/7/2011    0.85     1   2   1/10/2011   -0.14


title: Determine concrete important factors in pca

text: after using the R function pcromp () I get how many factors might be important. How to determine what are these factors? Which statistical tool could I use? I appreciate your reply!


title: How to build an "equilibrated histogram"?

text: I would like to create an "equilibrated histogram" with roughly the same number of data points in each bin. The second restriction I have is that I should have only 4 bins. Given the following list of numbers below, how can I achieve that? -2.153, -1.732, -1.699, -1.559, -1.355, -1.306, -1.151, -1.129, -0.636, 0.4085, 0.5408, 0.5731, 0.5842, 0.6206, 0.8175, 0.8274, 0.8710, 1.3214, 1.5552, 2.2342 Thanks!


title: Taking advantage of many pre-treatment measurements

text: I am planing a pre-post treatment-control design study with a large number of pre-treatment measurements. I have subjects divided into a control group and a treatment group. For both groups, I will collect hourly data for one year prior to the start of the treatment and then continue collecting data for another year. This will yield approximately 9000 pre-treatment measurements and 9000 post treatment measurements for each subject. The treatment is something that cannot be stopped once it is started, so a cross-over design could only be of the form AA/AB, which won't take advantage of the benefits of that type of design. The psychological and bio-medical literature suggests using an ANCOVA model, where the pre-treatment data is used as a covariate in the model. Putting 9000 covariates in a model seems totally ridiculous. Also, reducing the pre- treatment data to a summary statistic doesn't take advantage of the large number of measurements. I'm sure that this must have come up before, any ideas? References to published results would be especially helpful.


title: Comparing regression models

text: I recently fit 4 multiple regression models for the same predictor/response data. Two of the models I fit with Poisson regression.               model.pois <- glm(Response ~ P1 + P2 +...+ P5, family=poisson(), ...)     model.pois.inter <- glm(Response ~ (P1 + P2 +...+ P5)^2, family=poisson(), ...)      Two of the models I fit with negative binomial regression.               library(MASS)     model.nb <- glm.nb(Response ~ P1 + P2 +...+ P5, ...)     model.nb.inter <- glm.nb(Response ~ (P1 + P2 +...+ P5)^2, ...)      Is there a statistical test I can use to compare these models? I've been using the AIC as a measure of the fit, but AFAIK this doesn't represent an actual test.


title: Multi class LDA vs 2 class LDA

text: The problem of designing a multi-class classifier using **LDA** can be expressed as a **2 class problem** (one vs everything else) or a **multi- class** problem. Why is it that in certain cases **Multi-class LDA classifier** out-performs **2 class LDA** (one vs everything else) or **vice-versa**.


title: controlling for a highly correlated variable

text: In a paper looking at hand movements and time taken to complete a task it appears that the ratio of movements over time is constant (i.e., movement and time are highly correlated) yet they say in the paper they can control independently for both variable - how is this possible? > There is a strong relationship between the time taken and number of hand > movements made (Spearman coefficient 0.79, P <0.01). This has been > demonstrated with ICSAD before. Therefore, why not just time the procedure > with a stopwatch? This is answered when we apply partial correlation > coefficient tests. When controlling for time, the number of movements made > significantly compares with surgical experience and global score > (correlation coefficient −0.44 and 0.56, respectively, P <0.01 for both). > However, when controlling for movement, the time taken had no such > relationship with experience and global rating (correlation coefficient > −0.02, P = 0.9; 0.10, P = 0.8, respectively), suggesting that operative > speed is secondary to economy of hand movement. Datta V, Chang A, Mackay S, Darzi A. The relationship between motion analysis and surgical technical assessments. Am J Surg. 2002 Jul;184(1):70-3.


title: Split-half reliability on binary scores

text: I'm very new at statistics, so excuse me if my questions are to simple or obvious. I have a dataset with response times for two conditions. These response times were converted into a binary score of either 0 (absence of the effect) or 1 (presence of the effect). The experiment is designed in 7 blocks. For each block participants receive a 1 or 0. The total score for each subject can range from o to 7. So, I want to do a split-half reliability to see if the criteria I used to evaluate the presence of the effect is a reliable indicator of what I'm studying. At first I just split the column of 1s and 0s in half and correlated it, but that give me a nonsense value (something like 0.0009) so I'm missing something here. Conceptually I know I should correlate two parts of the test that are supposedly measuring the same thing, but I'm having difficulties at figuring out exactly how to do this for this dataset. Any help is much appreciated.


title: How to check whether a sample is representative across two dimensions simultaneously?

text: I'm attempting to develop a standardized method to check whether one set of locations are representative of a larger set. In this particular case, I'm attempting to look specifically at their geographical representativeness. One method is to look at a two sample t-test for latitude and longitude * independently*, but that clearly ignores the possibility that there could be correlation between the values. Another option is to look at a categorical grouping of location (eg state, market, or any other gridding of the network), and use a chi-squared test. Neither of these strike me as optimal, however. Is anyone familiar with a test that can check bias of a sample based on two dimensions simultaneously? Any thoughts would be greatly appreciated.


title: Discrimination between multivariate populations when variables are non-normal?

text: I'm analyzing chemical profiles of amphetamine samples, and each sample has 26 numerical variables describing the relative quantity of some chemical constituent. Only a few of the 26 variables have pretty Gaussian distributions, the rest are amorphous (e.g. one very large "0" bin and a small Gaussian-like hill to the right). I'd like to test whether populations differ by geographical origin (or some other ordinal variable), but I'm not sure how to test against this null hypothesis because the underlying assumption of multivariate normality is violated. Can I perform a Hotelling's T-square test via permutation, or am I looking for something like the Kullback-Leibler divergence? My training is in biology and most certainly did not prepare me for this, so I apologize in advance for having missed obvious things. Thanks a lot.


title: method of moments with variance=$\sigma^2$

text: I am trying to estimate the value of a parameter by equating variance from a distribution to the sample variance... i.e. using method of moments estimation. Would it better to use the variance formula with $1/(n-1)$ in the denominator or $1/n$? Why is one better than the other? thanks * * * I'll give an example: say I have been told Random variables from the 'geometric distribution' have a $\sum_{1}^n(X_i - \bar{X})^2$= 10. Then in order to estimate the parameter 'p' where $\frac{1-p}{p^2}$ is the variance of the geometric distribution. do I in this circumstance use: $\frac{10}{n-1}=\frac{1-p}{p^2}$ and then given 'n' solve quadratic for p. or should I be using $\frac{10}{n}$ instead on the LHS.


title: stochastic vs deterministic trend/seasonality in time series forecasting

text: I have moderate background in time series forecasting. I have looked at several forecasting books, and I don't see the following questions addressed in any of them. I have two questions:   1. How would I determine objectively (via statistical test) if a given time series has:     * Stochastic Seasonality or a Deterministic Seasonality     * Stochastic Trend or a Deterministic Trend   2. What would happen if i model my time series as a deterministic trend/seasonality when the series has a clearly stochastic component? Any help addressing these questions would be greatly appreciated. Example data for trend:               7,657     5,451     10,883     9,554     9,519     10,047     10,663     10,864     11,447     12,710     15,169     16,205     14,507     15,400     16,800     19,000     20,198     18,573     19,375     21,032     23,250     25,219     28,549     29,759     28,262     28,506     33,885     34,776     35,347     34,628     33,043     30,214     31,013     31,496     34,115     33,433     34,198     35,863     37,789     34,561     36,434     34,371     33,307     33,295     36,514     36,593     38,311     42,773     45,000     46,000     42,000     47,000     47,500     48,000     48,500     47,000     48,900


title: Heteroscedasticity in the context of general additive model

text: For my given data (identity link), my target is to check a political variable, i.e., I intend to check, if a certain impact, given by some treatment at a certain time, affects the intercept and slope of another dependent variable. To get a feeling of how the continuous, explanatory variables will impact the dependent one, I fitted a general additive model (GAM) with a b-spline smoother. Since taking the log of a independent variable, a parametric reformulation seems appealing for later analysis, calculation of the area beneath the curve, I logged that particular variable. As it now seems, plots of the squared residuals vs the dependent variable, heteroscedasticity seems to be an issue. Since I don't want to log the dependent variable too, though it would help a little, I wonder what other ways are there to deal with heteroscedasticity in the GAM context. Is it just like in the OLS context where one uses heteroscedastic and autocorellation robust standard errors?


title: Is this a repeated measures design or not?

text: How do you describe, or what do you call, a test that uses two factors as independent variables but uses a dependent variable that is a difference between two measures taken repeatedly from the same individuals? Usually, repeated measures means that multiple independent variables are measured on the same individuals. Here, it is the dependent, outcome variable that is the repeated measure. But, since the difference is taken for each individual, is this simply a two-way ANOVA?


title: GLM with two related predictors (X and X-squared)

text: I am running a general linear model in which I have two predictors: X and X-squared. I entered both these predictors in my analysis because I think X might explain the variance in the outcome measure partially linearly and partially in a quadratic fashion. Obviously, there is multicollinearity in this example. However, I was wondering if there are maybe some reasons why it is not a good idea to put both these predictors in my model.


title: Least squares minimization with unweighted data

text: A couple of questions regarding least squares minimization with unweighted data: 1) Assuming you do not have any estimates of errors/weights for data, is there any way to provide a goodness of fit? If the errors are known, then the $\chi^2$ can be used to do a goodness of fit test, assuming that the errors are gaussian. 2) What are your thoughts on weighting the data by the actual y values, i.e. $w_i=\frac{1}{y_i}$. This seems reasonable in my case, as the data should have larger errors for larger y values. Even so, without accurate estimates of the errors, the reduced $\chi^2$ is still not likely to be $\sim1.0$, so the goodness of fit tests don't mean anything. So basically, using unweighted data, we can find best fitting parameters, but thats about it?


title: In statistics what does NA stand for?

text: I understand that NA means data is missing, null or not present. But what do the letters NA stand for? "Not Available"?


title: How to treat a time rank in logistic regression?

text: For an analysis of user attrition I want to perform a logistic regression, using task difficulty and point in time (100 trials) for a math task. My guess is that high difficulty at a late point in time leads to increased probability to attrite, while high difficulty in early stages will not lead to a negative effect.   In my first try I simply included this time rank as a numerical predictor, but this obviously violates the normality assumption. But does it make sense to enter the 100 trials as order factors each? I also thought about binning these trials in order to reduce factor levels but wanted to ask for some opinions first.   I would be very grateful for any help with this problem.


title: Numeric methods to estimate curve parameter from 2 input data with Gaussian noise

text: I have:   * a function of the form $S=M f(t, \theta)$ where $\theta$ is a variable, $t$ and $M$ are parameters.   * two observations $(\theta_1, S_1)$ and $(\theta_2, S_2)$ in which $S_1$ and $S_2$ are observed values for $t$, and $\theta_1$ and $\theta_2$ are the observed values of $\theta$. Furthermore, we know that my $S_1$, $S_2$ contains noise, which are normal distribution $N(0, d_1)$ $N(0,d_2)$ with unknown variance $d_1$, $d_2$. Now I want to estimate $t$. Mathematically, if $S_1$ and $S_2$ are observed precisely, we have $$S_1= M f( t, \theta_1)$$ $$S_2= M f( t, \theta_2)$$ Then, by division we have $\frac{S_1}{S_2}=g(t,\theta_1,\theta_2)$ following which the parameter $t$ can be estimated by a Mathematica tool like Findroot w.r.t. the parameter $t$. However, the problem here is the noise of normal distribution. The goal is to find a $t$ as precise as possible, in the sense of some probability-theoretic characteristics. The preferred way is to use a Matlab or Mathematica function, although I don't think they have a direct solution. I am totally new to numeric analysis, and totally new to Matlab. My naive idea is that:   1. Should the first step be removing the parameter $M$?   2. If so, there are infinitely many ways to do that, among which we can list      * $\frac{S_1}{S_2} = \frac{f(t, \theta_1)}{f(t,\theta_2)}$     * $\frac{S_1+S_2}{S_1-S_2} = \frac{f(t, \theta_1)+f(t,\theta_2)}{f(t, \theta_1)-f(t,\theta_2)}$     * $ \frac{S_1^2 }{S_2^2} = \frac{f(t, \theta_1) ^2 }{f(t,\theta_2) ^2}$ It is clear that all these variations are equivalent mathematically if $S_1$, $S_2$ are observed precisely. However, they should yield different $t$ (because of the Gaussian noise in $S_1$ and $S_2$) My questions are   1. Do you think my first step should be to remove $M$ using the foregoing approach?   2. Is there some kind of combination of $S_1$, $S_2$ that gives better results than others?   3. What would be the best way to estimate $t$ of this problem? (Maybe in Matlab, or Mathematica) Thank you. I hope I was clear.


title: Using Bayesian model diagrams to present both model description and results (posteriors)?

text: The model diagrams in "Doing Bayesian Data Analysis", John Kruschke creates diagrams like this: ![enter image description here](http://i.stack.imgur.com/bGw94.jpg) To represent The following BUGS/JAGS code: ![enter image description here](http://i.stack.imgur.com/yw4vT.jpg) He discusses this representation in his related blog post, Graphical model diagrams in Doing Bayesian Data Analysis versus traditional convention I just wrote one out for my rather complicated model and had a real moment of clarity. The diagrams are generic representations of probability density functions- they don't really reflect the flat priors stated in the model (although technically, the x-axes are not labeled). I haven't actually read the book, but I think it would be even more useful to use this diagram to present both the model structure _and_ the results, i.e., replace these generic distributions with the posterior probability density of each parameter (e.g. $\textrm{N}(\textrm{figure of }M_0,\textrm{ figure of }T_0)$ (rather than a single distribution to represent both parameters as is the case below. I have three questions:   1. Are there any potential technical (statistical) issues with this approach? I only ask because I haven't seen results presented alongside the model description in this way.   2. Are there any other suggestions for how to make this do a great job of communicating results? One idea that I had was to present an indexed beta (e.g. a random effect for 1...n categorical treatments, in the model above it would be `beta1[i]`) as overlapping densities (one for each treatment effect)?    3. How could I make this presentation be more intuitive to readers not familiar with hierarchical modeling? (I think this is the most intuitive presentation that I have found so far, but there may be room for improvement)


When statistics is solving all the problems, why does one study the fuzzy theory? Is there any problem which can be solved only by fuzzy theory but not by statistics? Please kindly answer with examples.


title: Dealing with 0,1 values in a beta regression

text: I have some data in [0,1] which I would like to analyze with a beta regression. Of course something needs to be done to accommodate the 0,1 values. I dislike modifying data to fit a model. also I don't believe that zero and 1 inflation is a good idea because I believe in this case one should consider the 0's to be very small positive values (but I don't want to say exactly what value is appropriate. A reasonable choice I believe would be to pick small values like .001 and .999 and to fit the model using the cumulative dist for the beta. So for observations y_i the log likelihood LL_iwould be                if  y_i < .001   LL+=log(cumd_beta(.001))      else if y_i>.999  LL+=log(1.0-cum_beta(.999))      else LL+=log(beta_density(y_i))      What I like about this model is that if the beta regression model is valid this model is also valid, but it removes a bit of the sensitivity to the extreme values. However this seems to be such a natural approach that I wonder why I don't find any obvious references in the literature. So my question is instead of modifying the data, why not modify the model. Modifying the data biases the results (based on the assumption that the original model is valid), whereas modifying the model by binnning the extreme values does not bias the results. Maybe there is a problem I am overlooking?


title: List of likelihood-based classification techniques

text: This is a basic statistical pattern recognition question. I'm aware of **LDA classification** , **Naive Bayes Classification** techniques which give output as a **likelihood** (of data belonging to a certain class). What are other likelihood-based supervised classification techniques? Can all of the supervised classification techniques be expressed as **likelihood** or **likelihood ratio** in order to classify given data? Can anyone point me to articles of notes which can help me out?


title: Distance function for categories in K-means

text: How to define a distance function when euclidean distance doesn't apply? For instance, say I have some data involves nationality. I'll probably assign a number to each nation, but for nations that have smaller difference in numbers doesn't mean that they are more prone to be in the same cluster as nations that have bigger difference in numbers. Is it make sense if I just define a function that return 0 if two nations are the same, and return some positive integer otherwise? If so, how big that positive integer should be?


title: Optimal winsorizing cutoff for stratified proportional to size sampling

text: I'm trying to use Kokic & Bell 's ideas to find the cutoffs for the following population estimator $\sum_{h}\frac{1}{n_h}\sum_{i}\frac{M_{hi}}{z_{hi}}\frac{1}{m_{hi}}\sum_j y_{hij}^*,$ where we use stratification and proportional to size cluster selection and a subsampling, and $y_{hij}^*$ stands for the winsorized value. Is there any reference for this sampling design? It'd be helpfull if someone shares some experience in the problem of finding the cutoffs that minimize the mean squared error.


title: Decision Tree as variable selection for Logistic Regression

text: I have to do a Logistic Regression, and have to use a subset of the variables. I received this "tip": _do a Decision Tree first, and use the most relevant variables in the Logistic Regression_. Is this a valid technique? If not, what problems does it face? PS: This tip is also given here for a standard regression. **UPDATE** : First of all thanks for all the answers. I'm definitely going to try the LASSO model. But, the algorithm must be simple, and easy for the IT to implement. It's a "millions of lines" kind of database, and the model needs to recurrently be estimated. Does it change the approach?


title: What is the Fisher information for the truncated poisson distribution?

text: The zero-truncated poisson distribution has probability mass function: $$P(X=k) = \frac{e^{-\lambda}\lambda^k}{(1-e^{-\lambda})k!}$$, $k=1,2,...$ And the expectation of the truncated Poisson distribution via MLE is given as $\frac{\lambda}{(1-e^{-\lambda})}$ According to this document (pages 19-22) the Fisher Information is given by $$I(\theta) = \frac{n}{(1-e^{-\lambda})}\left[\frac{1}{\lambda}-\frac{e^{-\lambda}}{(1-e^{-\lambda})}\right]$$ How is this derived? =========================================================================== Here is my best attempt thus far (but still the wrong answer at the end): The likelihood is given by: $$L(\lambda) = \prod_{n=1}^n\frac{\lambda^{x_i}e^{-\lambda}}{x_i!(1-e^{-\lambda})}$$ and thus the log-likelihood is: $$l(\lambda)=-\ln(x_i!)+\ln(\lambda)\sum_{i=1}^nx_i-n\lambda-n\ln(1-e^{-\lambda})$$ Differentiating with respect to $\lambda$ we get the first derivative: $$l'(\lambda)=-n-n\ln(1-e^{-\lambda})+\frac{1}{\lambda}\sum_{i=1}^nx_i$$ and the second derivative with respect to $\lambda$ is: $$l^"(\lambda)=n\frac{e^{-\lambda}}{(1-e^{-\lambda})^2}-\frac{1}{\lambda^2}\sum_{i=1}^nx_i$$ The Fisher Information is given by $$I(\lambda)=E[-l^"(\lambda)|\lambda]=E\left[-\left(n\frac{e^{-\lambda}}{(1-e^{-\lambda})^2}-\frac{1}{\lambda^2}\sum_{i=1}^nx_i\right)|\lambda\right]$$ $$\Rightarrow -n\frac{e^{-\lambda}}{(1-e^{-\lambda})^2}+\frac{1}{\lambda^2}\sum_{i=1}^nE[x_i|\lambda]$$ Which yields my incorrect Fisher Information: $$I(\lambda)=-n\frac{e^{-\lambda}}{(1-e^{-\lambda})^2}+\frac{n\lambda}{\lambda^2}=n\left(\frac{1}{\lambda}-\frac{e^{-\lambda}}{(1-e^{-\lambda})^2}\right)$$ What have I done wrong with the expectation?


title: Derivation of Rayleigh-distributed random variable

text: I only have a uniform distribution function between [0,1]. And from this distribution, I should generate a sequence of Rayleigh distributed random variable using some software. Anyhow, I was able to finish the problem using the formula from the Wikipedia article: $$(1)\;\;\;\;X=\sigma\sqrt{-2\ln(U)}$$ However, there is only one thing I could not understand. I tried many times to derive formula $(1)$ using the Inverse transform sampling method, but I could not. Can someone show me the steps of how $(1)$ is found?


title: Randomly picking from $n$ choices roughly $n$ times. What's the resulting frequency distribution called?

text: Not sure of the best way of phrasing this question, but I'll give it a go. If I were to randomly choose whole numbers between 1 and $n$ a significant number of times relative to $n$ (say, $m$, where $m$ is something like 70% of $n$) and then looked at the distribution of frequencies of outcomes, I believe you would get a distribution that wasn't flat. I.e. a very small number of outcomes would come up >1 times, the majority would come up 1 time, and some outcomes would come up 0 times. Obviously, this would be affected by the relative size of $m$ to $n$: If I were to choose a only 10 numbers between 0 and 1,000,000 ($m$=10, $n$=1,000,000), I would expect the vast majority of numbers to have a frequency of 0, with the rest 1. I have questions:   1. Is this correct?    2. How would you calculate the expected frequency of the most frequent chosen number?   3. If you were to order the outcomes by frequency and then bar graph the frequency, does the resulting curve have a name, or any interesting properties? I hope this makes some sort of sense. If not, let me know and I'll try and think of a better way of explaining it. This all came from an original question which went something like this: "If I have a lottery machine with 15 million different lottery number outcome combinations, given that I've already made 10 million draws, what is the probability that the 10,000,001th draw has already been drawn?" **Edit:** I ran some monte-carlo simulations for $2<n<650$, with $m$=$n$ and the total number of experiments for each $n$ run 10,000 times and the results averaged. I logged the max frequency, the number of times a number came out 1 time, and the number of times a number comes out 0 times. Two interesting things: The number of numbers that come up 1 time closely matches the calculated $n p(k)$, and this is also close to the number of times a number comes out 0 times. I'm not sure it's obvious that these two should be the same, as $p(k)$ doesn't seem to make sense for $k=0$ (Can you do $m$ choose $0$?). Also, interestingly, a plot of $m$ against the average maximum frequency gives a curve that looks like this: ![](http://i.imgur.com/h0Ttp.png) Obviously some sort of log relationship. Next step is to vary $m$ as a percentage of $n$.


title: Bias in sampling for set intersections

text: Say I have 2 sets, $A$ and $B$ with $n_{A}$ and $n_{B}$ elements respectively, which I assume is known. I would like to estimate $| A \bigcup B |$ using samples of $\tilde{A} \subset A$ and $ \tilde{B} \subset B$. That is if $\tilde{A}$'s elements are uniformly sampled from $A$, and likewise for $\tilde{B}$, will $ \frac{| \tilde{A} \bigcup \tilde{B} |}{| \tilde{A} | + | \tilde{B} |}$ be an unbiased estimate for $\frac{| A \bigcup B |}{| A |+|B |}$? If not, is there some another estimator that will allow me estimate $| A \bigcup B |$ without bias?


title: Exponential moving average with sub-interval relevance / varying timeframe

text: I need to calculate an exponential moving average for a series of data. The intended sampling interval is fixed (say 1s) but the data stream has varying intervals (data intervals vary from 0.01s to 10s or so). The data is somewhat noisy (a random data sample would virtually never be on the average). My impression is that I cannot thus just take the most recent data sample at each interval, as that could easily lead to a misleading stat. I reckon I can somehow just calculate the average for each period and take that as the sample, but I'm not positive. Is there a standard algorithm that manages an exponential moving average on a time-variable data stream? * * * _I need to program this for a real-time system where it won't be possible to store a sample history. Nonetheless, I'm sure I can adapt any non-streaming algorithm to the streaming form._


title: Can SD be calculated from binned data? How to test non-random groups that are clearly different?

text: I DON"T have the underlying sample data (but could possibly request it). We are talking millions of data points (utility customers). It is actually the full population (not a sample) combed for those customers that were not in one rate class for the entire 2011 year (about 75% pass this threshold). Since the full data set is highly cumbersome, I received it in highly bracketed form (e.g. x,xxx customers used between 6,000 and 6,100 kWhs for the year, x,xxx customers used between 6,100 and 6,200 kwh). **1\. Can an approximation of standard deviation be calculated using this bracketed data? Is there a convenient way to do this in Excel? Fortunately, the brackets are narrow (there may be as many as a thousand).** I have two utility populations that are not at all like a randomly assigned control and test group for a medicine test. They are customers who qualify and sign up for a subsidy low income rate and those who don't. They are by definition, clearly different. The smallest usage difference I experienced was 2%. Despite what might test as fairly high variance of usage (some poor use a lot and some rich are frugal), given a test of 3 million customers, I predict the difference would test as highly significant. **2\. Is Welch's t-stat the appropriate test for these two populations?**


title: Preferred methods for graphing time-series data to present "averages"?

text: I have data from a load test of a web site with several thousand data points spread out over roughly 30 minutes (the values are the response time of the site in milliseconds). The values are spread out among the 30 minute range, but not at a constant rate (i.e. there may be a few milliseconds between some points, other points maybe at the same timestamp, etc). I'd like to present this data visually and chart it, but I'm not sure of the best method for doing so - there is a good amount of variance around any sort of concept of average values or a trend line. Is there any generally accepted best practices for methods on how to graph data of this type? I'm concerned about choosing a poor method for averaging/smoothing out the data and misrepresenting the data - such as underweighting some outlier values. I've played around with a line chart with the timestamps on the x-axis and the average of samples in the same minute on the y-axis. I'd also like to consider graphing a moving average of the data, but I'm unsure if I should be averaging datapoints in the same N minutes or a window of the last N points. I'd like to make sure that whatever choice I make would appear to be a rigorous representation of the data and not too amateur-ish. **Update:** below is a sample of what I have produced so far, each point on the chart is taken as the mean/median of all of the samples within the same minute (i.e. within 11:12:00.000 and 11:12:59.999). I included the number of samples per minute as a bar chart in the second half of the image to be able to show if any single points in the line chart look as outliers due to a small amount of samples, although aesthetically speaking I think the bar chart takes up way too much real estate for the amount of information it gives. ![](http://i.stack.imgur.com/4t1sT.png)


title: pgmm fitted.value

text: I'm running a dynamic panel model with lagged dependent variable, but seems after fitting the model, the fitted value not corresponds `log(RPK)` estimated:               d=pgmm(log(RPK) ~ lag(log(RPK),1) + log(GDPPR)+log(CPR)+log(NP)|lag(log(RPK),2:99),             data = no2, effect = "individual", model = "twosteps",collapse=TRUE,index=c('Flow','Year'))      fit=fitted(d)       In fact, the fit I obtained is much smaller than my initial `log(RPK)`, so I wonder maybe when using `fitted()`, it actually gives us some kind of difference? Is there anyone who could help me out about that?


title: Quick Question - Approximate distribution for sample mean?

text: ![enter image description here](http://i.stack.imgur.com/Ek05F.png)![enter image description here](http://i.stack.imgur.com/EOz4o.png) I am having issues answering part two. I think it is about CLT, correct me if I'm wrong. But how do you compute the distribution details from this ? Please help and thank you for all your contributions. Please answer it using the main question AND part one which is about the expression.


title: Good resources about reversible-jump MCMC sampler in R

text: Please help me to find out simple (I mean useful with little knowledge on Bayesian statistics) resources:   1. To learn MCMC and perform general linear or mixed models    2. To learn different sampling techniques including reverse jump sampler. Also, I want implement these methods using R base or R packages. I can use SAS and any other specially designed software, however I have R as preference.    3. Example codes to perform MCMC and reversible jump sampler.


title: Edge Probabilty Using KL-Divergence Code in Python

text: Its a little complicated question, so please bear with me. I am doing Image Segmentation using Swendsen-Wang method for Image Analysis./ ( _stat.fsu.edu/~abarbu/papers/jcgs.pdf_ ) I have to calculate the **edge probability** as a good similarity measure between the intensity models of the atomic regions. The paper states that(first of all) - > As an approximate model for each atomic region v ∈ V , we choose a 15-bin > intensity histogram h normalized to 1. So, I have created the histogram of an image, using this python code - > k = cv2.imread("vpython/the.png",0) > > hk = cv2.calcHist([k],[0],None,[15],[0,256]) Further, the paper states that we have to choose the edge probability as - > q(i,j) = exp{−1/2 (KL(h(i) || h(j)) + KL(h(j) || h(i)))} where KL() is the Kullback-Leibler divergence between the two histograms. I have tried coding it as - > for i in range (int(k.size)): > > for j in range (int(k.size)): > > q(i,j) = q(i,j) + math.exp(((hk(j)-hk(i))*math.log(hk(i)/hk(j)))/2) which is obviously not correct, but that is what is my understanding of this problem. I am getting this error - > File "hist1.py", line 14 > > q(i,j) = q(i,j) + math.exp(((hk(j)-hk(i))*math.log(hk(i)/hk(j)))/2) > > SyntaxError: can't assign to function call Can anyone help me with this? Thanks in advance!


title: Is Naive Bayes fine for simple "Suggested products" solution?

text: I have the opportunity to develop a simple solution to display suggested products on an e-commerce website. My first choice was to use Naive Bayes due to its simplicity of implementation. However, I got to this state: $P(buys|product.type,product.brand) = \frac{P(buys) * P(product.type|buys) * P(product.brand|buys)}{P(product.type,product.brand)}$ As we compare P(buys) and P(does not buy), the denominator is a constant and we can compare P(buys) and P(does not buy) based on their numerator only. Determining the probabilities $P(product.type|buys)$ is easy by looking at what the user bought in the past. However, I'm not sure how I would determine $P(buys)$. Is it the number of products bought divided by the total products available in the shop? But the number of products available varies with time. Is it fine to take the current number of products, which will be different from the number of products that were available when the user made his last order (an extreme example: the user can have ordered 3 of the 5 products that were available at launch, but now there are 10000 products).


title: Controlling overfitting with random forests for very high dimensional data

text: I'm using the randomForest package in R to analyse a genetic dataset of ca 100 samples with 10.000 genes. The samples are grouped into 5 classes, the smallest being only of 5 samples. I'm ultimately interested in the proximity matrix, confusion matrix and the variable importances, which is why I don't wish to start by reducing the dataset with PCA or similar. Rather, I'd tune the model using OOB estimates to handle the high dimensionality and variable class - if possible at all! What's the best parameters to tune for such high dimensional data? After some initial tests, raising the minimum end-node size seems to improve OOB scores. I was thinking the main parameters of interest should be:   1. Number of trees in forest (ntree)   2. Number of features sampled (mtry)   3. End-node size (nodesize) I don't know the best way to control for very small class sizes. The small classes are never predicted correctly in my initial tests. What's the best way to solve this problem? Simply to run all possible combinations of parameters, and pick the best based on OOB estimates, or is there a smarter way?


title: Why linear transformation can improve classification accuracy when the dimensionality of data is high?

text: Let $X$ be an $m\times n$ ($m$: number of records, and $n$: number of attributes) dataset. When the number of attributes $n$ is large and the dataset $X$ is noisy, classification gets more complicated and the classification accuracy decreases. One way to over come this problem is to use linear transformation, i.e., perform classification on $Y=XR$, where $R$ is an $n\times p$ matrix, and $p<=n$. I was wondering how linear transformation simplifies classification? and why classification accuracy increases if we do classification on the transformed data $Y$ when $X$ is noisy?


title: Probability of Event from a table

text: I have the following table: ![SPSS output](http://i.stack.imgur.com/sWfCn.png) From here I want to answer the following question: P(CHILD CAN RECOG BEG SOUNDS OF WORDS – YES | CHILD RECOGNIZES LETTERS - Most of them). From my understanding of this I am saying 81.1%, am I reading this correctly? Child Can Recog (Recognize) Beg (Begining) Sounds of Words. So what I am trying to find from the data, is: The Probability that a Child can recognize beginning sounds of words (1 - Yes) given that the child recognizes most of the letters. So what I did was go to the Row 1-Yes, then the row in that area labeled % within PE2 - Child recognizes letters and then over to column 2 called 2 Most of Them and picked the 81.1% **UPDATE** I ended up doing the following:                                     RECOG LETTER - MOST | RECOG LETTER OTHER | TOTAL     RECOGNIZE WORD YES          361           |       858          | 1219     RECOGNIZE WORD NO           84            |       232          | 316     TOTAL                       445           |       1090         | 1535          361/1535 = 0.24, 858/1535 = 0.56, 84/1535 = 0.05, 232/1535 = 0.15      Which gave me the following:                                      RECOG LETTER - MOST | RECOG LETTER OTHER | TOTAL     RECOGNIZE WORD YES           0.24          |       0.56         | 0.8     RECOGNIZE WORD NO            0.05          |       0.15         | 0.2     TOTAL                        0.29          |       0.71         | 1          P(RECOG LETTER - MOST) = 0.29          P(RECOGNIZE WORD YES) = 0.8          P(RECOG LETTER MOST AND RECOG WORD YES) = 0.24          P(RECOG LETTER MOST OR RECOG WORD YES) = 0.24 + 0.05 + 0.56 = 0.85      So then I did the following:               P(RECOG WORD YES | RECOG LETTER MOST)     = P(WORD YES AND MOST) / P(WORD YES)     = 0.24 / 0.8     = 0.3


title: nonlinearTest() in S-plus vs R

text: Is there any similar function in R does the same thing as nonlinearTest() in S-plus? I am currently working on determine the number of regimes,delayed parameter d and threshold values for the threshold value. Or otherwise where can I download a student S-plus? Do I have to use the university computer in order to free download?


title: Treatment of missing values introduced by padding lagged variables

text: In the case of a linear regression with lagged independent variables, what are the techniques for dealing with the NA values introduced by padding lagged variables (since t < 0 values do not exist)? I ask because I'm implementing a random forest with lags and have to decide what to do about these missing values. Replacing them with 0 seems misguided, but imputing backwards doesn't make much sense either. I suppose I could insert a time mean, or an overall mean for the variable if I have panel data... Edit 1: Looking around more, SAS lets you choose between replacing NA with a time mean, overall mean, or zero: http://support.sas.com/documentation/cdl/en/etsug/60372/HTML/default/viewer.htm#etsug_panel_sect015.htm Any idea on the costs and benefits of these different approaches?


title: How to choose training sets for textual data?

text: I have textual data scraped from a classified ads site and I'm interested in classifying the posts into two categories: "X" and "not X". I'm stuck with coming up with a good training set and feature selection. I assumed that a good start for a "corpus" of text about X would be its product page from the manufacturer (after all, a seller would have to use _some_ technical terms) and pertinent technical jargon as the "features". I've used NLP to process the text to get rid of common/filler/free-style text,leaving me with just 3-5 keywords per target post. Does this approach seem right? Also, what would I use as training for "not X"?


title: Variance pooling when sample size is a predictor

text: Suppose that I am building a hierarchical model of performance and the data is hierarchically structured (e.g., multiple customers rating a single salesperson). I might want to use variance pooling in this situation to allow individuals who have not had many customers yet borrow information from those who have had more customers. Then again, those with fewer customers may have less experience, and may generate fewer leads because they are poor quality salespeople in general. Here, the sample size is correlated with one or more possible predictors. How to account for this in hierarchical models, especially Bayesian ones? Thanks.


title: linear model, given y, calculate the confidence interval for x

text: I fit the model `y = a + b * x`. And 95% CI for estimates of a and b are (a1, a2), (b1, b2), respectively. If we have a new observation `x0`, then the estimated response is `a_hat + b_hat * x0`. And the 95% CI is `(a1 + b1 * x0, a2 + b2 * x0)` Case 1, now given `y = 0`, I want to estimate the responding `x`. I guess the confidence interval for `x` is `(a1 / b2, a2 / b1)`. Am I right? Case 2, if we fit a model `y = a + b * x + c * z + d * x * z` , `a_hat<0, b_hat,c_hat,d_hat>0`, Given y = 0, how to calculate the confidence interval for `x + z`?


title: Which Xcorr to choose in a concrete situation: normalized/not, biased/unbiased?

text: I have 2 data sets, one is binned and another is not, 9000 values each. Something like this: `a = [23 45 88 9 0 29 30 60 0 2 10 100 0 60 90] b = [0 0 0 1 1 0 0 0 0 0 0 1 0 0 1]` Actually, the `b` is the set of events (yes/no) which occurs during a continuous change of the intensity of a "signal" (the signal has values [0–100]). I performed the cross-correlation, because I want to find if the occurrence of the events has something in common with the changes of a signal. So I did:               maxlags = 300; [c, lags] = xcorr(a-mean(a), b-mean(b), maxlags, 'coeff')      And I also created a 10-times shuffled null-hypothesis:                   urge_shuffled=urge_new;     a_sh = a;     b_sh=b;     all_corr_rand_norm=[];     for k=1:10         a_sh1 = a_shuff(randperm(length(a_sh)));         a_sh = a_sh11;         b_sh1 =b_sh(randperm(length(b_sh)));         b_sh=b_sh1;         [c_sh,lags]=xcorr(a_sh-mean(a_sh),b_sh-mean(b_sh),maxlags, 'coeff');         all_corr_rand_norm (k,:) =c_sh;     end         corr_rand=mean(all_corr_rand_norm);`      and also calculated the _t_ -score values. So finally I've got such a thing: ![](http://s30.postimg.org/iywttvlfl/Bildschirmfoto_2014_06_30_um_19_00_54.png) At first I tried to do the correlation just with `xcorr(a,b, maxlags, coeff)` without subtracting the mean, but in this case it gives me a very strange null hypothesis (it's in red; the real correlation is in blue; the changes are more evident with bigger maxlags): ![](http://s24.postimg.org/n2an9ah0l/Bildschirmfoto_2014_06_30_um_18_56_31.png) It's clearly a triangle. I did some research on the internet and even found some similar questions, but the answers were that this is because of "zero- padding", and that's why we should subtract the mean before the correlation. Here is my 1st question: **why in some cases do people get a nice null- hypothesis even with`xcorr(a,b)` – what does it depend on?** The second question is **whether or not it's right in my case to subtract the mean** – basically, it already corresponds to the `xcov` function of MATLAB, and this is already the covariance... (I also have a colleague who suggested I use the `xcorr` and afterwards to subtract the mean:               [c, lags]= xcorr(a,b, maxlags, 'coeff')     x_crr=c-mean(c)      —but I don't really see why I should do it, and besides, it still gives me the same strange triangle-shaped null-hypothesis. The 3rd question is: **are the cross-covariance and cross-correlation so different in my case?** Yes, I've already read all that I could find about them both, but the phrases, "The covariance shows how 2 functions variate together when the cross-correlation shows their correlation," didn't help too much. If something varies with a specific lag before/after the other function, isn't it the correlation? Also, while searching for the explanation of my triangulated null-hypothesis, I saw somebody's advice to use the "'unbiased' cross-correlation", while others strongly suggested to not do this. So the 4th question is, **when do we actually use this unbiased/biased`xcorr`?** I've found again mostly "great" descriptions that the "biased cross-correlation is biased, while unbiased is not." There were some better answers, of course, but I still don't get it. And, the last question: in my case, as I have a significant _t_ but a weak correlation between both sets of values, can we affirm that the fluctuations of both sets of data are almost independent? I perfectly understand that it's a mix of common and probably easy questions about the cross-correlation, but the explanations that I’ve found on the internet usually sound dubious or are not really like explanations, but just some fast assertions (like "the unbiased means is unbiased") and I feel unsure about interpreting my findings.


title: Sufficient number of sample to learn Bayesian network?

text: I want to construct Bayesian network for a 800 genes(genes are my node/variables). I have only 30 cancer samples and 30 normal sample.so I want to create network for cancer samples and for the normal samples whether my data is reasonable to learn Bayesian network ?


title: Comparison of CPH, accelerated failure time model or neural networks for survival analysis

text: I am new to survival analysis and I've recently learned that there are different ways to do it given a certain goal. I am interested in actual implementation and appropriateness of these methods. I was presented with the traditional _Cox Proportional-Hazards_ , _Accelerated failure time models_ and _neural networks_ (multilayer perceptron) as methods to get survival of a patient given their time, status and other medical data. The study is said to be determined in five years and the goal is to give survival risks each year for new records to be given. I found two instances where other methods where chosen over the Cox PH:   1. I found " _How to get predictions in terms of survival time from a Cox PH model_" and it was mentioned that: > If you are particularly interested in obtaining estimates of the probability > of survival at particular time points, I would point you towards > **parametric survival models (aka accelerated failure time models)**. These > are implemented in the survival package for R, and will give you parametric > survival time distributions, wherein you can simply plug in the time you are > interested in and get back a survival probability. I went to the recommended site and found one in the `survival` package \- the function `survreg`.   2. Neural networks were suggested in this comment: > ... One advantage of neural net approaches to survival analysis is that they > do not rely on the assumptions that underlie Cox analysis... Another person with the question " _R neural network model with target vector as output containing survival predictions_" gave an exhaustive way of determining survival in both neural networks and Cox PH. The R code for getting the survival would be like this:              mymodel <- neuralnet(T1+T2+T3+T4+T5~covar1+covar2+covar3+..., data=mydata, hidden=1)     compute(mymodel,data=mydata)        3. I went to the R forums and found this answer in the question " _predict.coxph and predict.survreg_ ": > Indeed, from the `predict()` function of the `coxph` you cannot get directly > "time" predictions, but only linear and exponential risk scores. This is > because, in order to get the time, a baseline hazard has to be computed and > it is not straightforward since it is implicit in the Cox model. I was wondering if which of the three (or two considering the arguments over Cox PH) is best for getting survival percentages for time periods of interest? I am confused which of them to use in survival analysis.