title: How to cure Poor Fit Indices for CFA?

text: I am facing problem in getting the good fit indices. I have only 163 respondents with a total of 4 variables and 70 items ( **A:** 5 components, 3-4 items each; **B:** 3 components, 8 items each; **C:** 4 components, 4-5 items each, **D:** 11 items). The initial values for fitness indices were: ChiSq P-value = .000; RMSEA = .071; GFI = .586; SGFI = .558; CFI = .675; TLI = .663; NFI = .488; ChiSq/df = 1.815 Then, my supervisors suggested me to run EFA for the model but it ended up the items in the same group were being distributed to different groups, where the grouping of those items doesn't make sense at all as I adopted these constructs from past literature. The number of factors was not the same as the one proposed in past literature. **Q1: Should I modify the model according to the EFA result?** Without modifying the model, I tried to do CFA for each construct separately and deleted items by referring to Modification Indices (MI). And I ran CFA for entire model again by deleting what I've deleted in the last step. However, the fitness indices are still poor. If I further delete items with highest MI, the result is even worse. Then, I started over again by using Cronbach's Alpha to identified the items with lower reliability and deleted it. And I repeated the same steps by running CFA for individual constructs, deleting items according to MI, running CFA again for whole model, and deleting items according to MI. I've deleted 17 items from the model. The final result is as below: ChiSq P-value = .000; RMSEA = .060; GFI = .690; SGFI = .660; CFI = .831; TLI = .822; NFI = .649; ChiSq/df = 1.587 **Q2: May I know is there any other way that I could do to improve the fitness indices? If the problem is with my data, would I have to give up using SEM for analysis?** Thank you!


title: Trying to figure out a formula and answer to probability question

text: If I go to a certain establishment 12 times in one year (let's just say once a month), and there's a 6.25% chance each time that I go in there that I would run into a person I know, what is the probability that I would run into a person I know at least one time during the course of those 12 occurrences? Can you tell me the probability and exactly how you derived that answer?


title: Log-uniform distributions

text: I am having some difficulty understanding what log uniform distributions are. Suppose that $\log X$ is uniformly distributed on the interval $[1,e]$. How do I describe $P(X=x)$? It seems like there is more probability mass on the lower numbers so that X itself is not uniformly distributed, but I am having difficulty formalizing this argument.


title: Forecast accuracy calculation

text: We are using STL (R implementation) for forecasting time series data. Every day we run daily forecasts. We would like to compare forecast values with real values and identify average deviation. For example, we ran forecast for tomorrow and got forecast points, we would like to compare these forecasts points with real data we will obtain tomorrow. I am aware that forecasts values and real data may not match most of the times, that is one of the reason we would like to keep track of how much accurate we are each day. Now we are trying identify what is the best approach to solve this problem? any help pointers would be appreciated. I looked at Measuring forecast accuracy question, but it seems it is related to comparing models rather than calculating accuracy with real values. I looked at accuracy function implementation in R, but confused with two questions: 1) Will it work on real data vs forecast data, because most of the tutorial saying as "test data" vs "forecast data" 2) It seems out put of accuracy function is array of values rather than % of deviation.


title: How do you pronounce "LASSO"?

text: Some people pronounce it "LAS-so" and some "las-SOO" (from the rope cowboys throw to catch cows). Which way is "the right one"?


title: full batch v.s. online learning v.s. mini batch

text: This is a question from a coursera course: Suppose we have a set of examples and Brian comes in and duplicates every example, then randomly reorders the examples. We now have twice as many examples, but no more information about the problem than we had before. If we do not remove the duplicate entries, which one of the following methods will not be affected by this change, in terms of the computer time (time in seconds, for example) it takes to come close to convergence? a) full-batch learning   b) online-learning where for every iteration we randomly pick a training case   c) mini-batch learning where for every iteration we randomly pick 100 training cases The answer is b. But, I wonder why c is wrong. Isn't online-learning a special case of mini-batch where each iteration contains only a single training case? Thanks,


title: How to get P-Value when t value is less than 1?

text: I am doing an assignment for class and I have hit a wall. I'm not sure what to do next. Using a sample with 12 discrete elements, I am trying to test a hypothesis that the population mean is 243 with significance level .05. Here is what I have so far (this has been verified): n = 12 x bar = 244.33 s = 12.382 t = (244.33 - 243) / (12.382 / sqrt(12)) = .3720 The next step should be to calculate the P-value, which I do by looking up .3720 on the t-table with df = n-1 = 11, however there are no values on the t-table less than 1. Therefore I don't know how to look up .3720. I feel like there is some crucial bit of knowledge I'm missing that is stopping me from completing this problem. I would really appreciate some assistance in knowing how to calculate the P-value from here. Thank you very much. UPDATE: I appreciate everyone's comments but I feel as though I'm getting a lot of complicated answers. This is a question out of my textbook and isn't supposed to require any special software or anything like that, just a simple calculator and a t-table. Maybe it would be easier to ask the question from the text ad verbatim: *A random sample of n = 12 frozen dinners of a certain type was selected from production during a particular period, and the calorie content of each one was determined. Here are the resulting observations. The stated calorie content is 243. Carry out a formal test of the hypothesis that the population mean = 243, with a significance level of .05. Observations:   255 244 239 242 265 245 259 248 225 226 251 233 *


title: What is the difference between forecasting based on ARIMA and logistic curve? R

text: I'm making a project connected with identifying the dynamics of sales. My database concerns 26 weeks (so equally in 26 time-series observations) after launching the product. This is what my database looks like: https://imageshack.com/i/0yyh6ij I want to make forecast based on S-curve for clusters of time-series. The main aim was to compare two methods of forecasting:   1. based on parameters of logistic curve   2. based on ARIMA However, I do not know how to compare these two methods = measure their performance. That's a plot with prediction based on S-curve http://imageshack.com/a/img850/6600/rzkp.jpg So my questions are:   1. How to measure performance=forecast errors based on logistic curve?   2. How to compare forecasting based on logistic curve and ARIMA - what is the main difference between these two approaches if I base on one variable - units_sold_that_week? I would be grateful for any explanation.


title: Interpreting residual deviance

text: I have never understood residual deviance, other than the fact that it is a number that's useful for calculating $R^2$ of a linear Gaussian regression model. Let's say R outputs a residual deviance of 2,000, and a null deviance of 6,000.... So $R^2$ is then $1 - \frac{2000}{6000} = 0.6$, and we can say that "60% of the deviance can be explained with the regression...." But is there a way to explain, in plain English, what the values "2,000" and "6,000" actually mean? E.g. "The distance from the mean is, on average, 2,000 units," or something (I know that's not right).


title: Difference between "common index" and "no common index"

text: I have been looking at the G*Power tool for Statistical Power Analysis but I am confused by the difference between two of the z-test options.               "Correlations: Two dependent Pearson r's (common index)"      and               "Correlations: Two dependent Pearson r's (no common index)"      Is this in relation to the test data? e.g. A group of 10 men and a group of 10 women were asked to rate a single index of 10 images (common index) A group of 10 men and a group of 10 women were asked to rate two separate indexes of images (no common index) Any help is much appreciated


title: Significance of coefficients in linear regression: significant t-test vs non-significant F-statistic

text: I'm fitting a multiple linear regression model between 4 categorical variables (with 4 levels each) and a numerical output. My dataset has 43 observations. R gives me the following $p$-values from the $t$-test for every slope coefficient: $.15, .67, .27, .02$. Thus, the coefficient for the 4th predictor is significant at $\alpha = .05$ confidence level. On the other hand, R gives me a $p$-value from an $F$-test of the null hypothesis that all my slope coefficients are equal to zero. For my dataset, this $p$-value is $.11$. My question: how should I interpret these results? Which $p$-value should I use and why? Is the coefficient for the 4th variable significantly different from $0$ at the $\alpha = .05$ confidence level? I've seen a related question, $F$ and $t$ statistics in a regression, but there was an opposite situation: high $t$-test $p$-values and low $F$-test $p$-value. Honestly, I don't quite understand why we would need an $F$-test in addition to a $t$-test to see if linear regression coefficients are significantly different from zero.


title: Interpret significant predictor in non-significant regression?

text: > **Possible Duplicate:**   >  How can a regression be significant yet all predictors be non-significant? In a simple linear regression with multiple predictors, is it valid to interpret the significance of individual predictors even though the overall model is not significant?


title: How can you have a non-significant multiple regression model w/ significant predictors?

text: > **Possible Duplicate:**   >  Not-significant F but a significant coefficient in multiple linear > regression   >  How can a regression be significant yet all predictors be non-significant?   >  Significance of coefficients in linear regression: significant t-test vs > non-significant F-statistic If in a multiple linear regression (enter method) the general model isn't significant (F>.05) but one of the predictors is significant (β<.05), should I consider it as a significant result?


title: Correlation is significant but regression is not- updated

text: I have 1 dependent variable (a type of behaviour) and 5 independent variables. The independent variables all have statistically significant positive correlations with the DV (p<.01). However, I conducted a regression analysis and only one of the independent variables is now statistically significant. After reading some other posts, I am assuming that this is due to multicollinearity between my IVs. My question is: how do I report this in my study? Is it ok to say there is a relationship due solely on the siginificant correlation? and include an argument on multicollinearity as a limitation to my regression analysis?


title: X and Y are not correlated, but X is significant predictor of Y in multiple regression. What does it mean?

text: > **Possible Duplicate:**   >  How can a regression be significant yet all predictors be non-significant? X and Y are not correlated (-.01); however, when I place X in a multiple regression predicting Y, alongside three other (related) variables, X and two other variables are significant predictors of Y. Note that the two other variables are significantly correlated with Y. How should I interpret these findings? X predicts unique variance in Y, but since these are not correlated, it is somehow difficult to interpret. I know of opposite cases (i.e., two variables are correlated but regression is not significant) and those are relatively simpler to understand from a theoretical and statistical perspective.


title: How can a regression be significant yet all predictors be non-significant?

text: My multiple regression analysis model has a statistically significant F value however all beta values are statistically non-significant. All the regression assumptions are met. No multicollinearity was found. Correlations among all predictors are all less than 0.60. What else might be the cause of the insignificant predictors?


title: Variables correlated but not significant predictors of dependent variable in multiple regression

text: Two of my independent variables are correlated with the dependent variable but in multiple regression they turned out to be insignificant predictors of the same. Is it possible? and how to report this difference in data interpretation?


title: Non-significant MANOVA with significant planned contrasts

text: > **Possible Duplicate:**   >  Not-significant F but a significant coefficient in multiple linear > regression   >  How can you have a non-significant multiple regression model w/ significant > predictors? My hypothesis was that social engagement would be associated with cognition. I set up a one-way MANOVA in SPSS to look at the association between social engagement (IV, 9 levels) and performance on 6 different cognitive tests. Because the IV had so many levels, I also set up a series of planned contrasts (equivalent to Helmert contrasts) so that I could see where the real difference lay. The overall ANOVAs came back as non-significant for all DVs bar one. However, for 3 of the tests, the same planned contrast came up as significant (p<0.01). How could the planned contrasts come up as significant if the overall ANOVA was not? And how do I report this? Please help!


title: how to interpret the result of significant regression of binary logistic yet all predictor is insignificant

text: I am performing research on Product attributes of mineral water toward purchase decision. However, the result of -2log likelihood test showed that my model is fit enough yet my result of Hosmer and Lemeshow is only 0.121. However, my Wald statistics showed none of my independent variables were statistically significant. How should I interpret this result?


title: Is it still meaningfull to watch the coeffiecent(s) of a regression equation even when it/they is/are not significant?

text: I am reading the regression analysis resutls in a paper which showed most of the coefficients b2 and b3 are not significant but the R2 is significant. Eventhough there aren't too many variable in this model(see the model below), that happend. How come? Is it still necessary to watch b2 and b3 and make conclusion like from b2 we can tell that there is a possitiv correlation between X and Y? The author did that but I am not sure whether one should do that if the coefficient is not significant at all. ![enter image description here](http://i.stack.imgur.com/PGCfb.png) ![enter image description here](http://i.stack.imgur.com/whvWK.png)


title: multiple linear regression model does not reach significance, however, there is a siginificant coefficient. SO is it reliable?

text: I have 30 observations, and I inculde 5 variables as independents using ENTER methods to conduct multiple linear regression. The results showed the multiple linear regression model does not reach significance (Not-significant F), however, there is a significant coefficient. In this situation, is the significant coefficient still reliable ? and, how make this situation happen?


title: How to explain significant correlation between independent and dependent variable but non significant regression test in the presence of mediation?

text: > **Possible Duplicate:**   >  Not-significant F but a significant coefficient in multiple linear > regression In my research, I use Optimism as the IV, Job Satisfaction as the DV and Work- Family Enrichment as Mediator. Path Model Analysis using AMOS shows that IV and DV are significantly correlated, but regression test was non-significant. What conclusion can we draw from this?


title: Two independent variables both correlate with the dependent variable, but none are significant in a regression analysis

text: I recently encountered a problem in relation to my study of conflicts and performance. The two types of conflict in the study, task and relationship conflict, both correlated with performance $(-.39, P < .01)$ and $(-.37, P < .01)$. When I put these in a regression, the regression model came out significant, and it explained $22.1\%$ of the variance $(F\ (3, 82) = 7.76,\ P < .01)$. None of the two IVs were, however, significant coefficients! There were no other variables in the regression analysis. How could this be?


title: How can a regression be significant yet both predictors be non-significant?

text: I know this question has been asked in a slightly different form here. But my question differs and because of the forum rules I can't post on that thread. I have 2 independent variables, n=32, highly significant multiple regression, but only the intercept/constant has a significant p value. Also, I've tested for multicollinearity, it's not in the readout from R below, but the variance inflation factor for the 2 dependent variables is 2.02. Durbin-Watson, added variable plot, P-P plot, standardized residuals plot, and cook's distance are all normal. Here's the R readout - please let me know what you think. ![enter image description here](http://i.stack.imgur.com/Pv6p9.jpg) SO: This topic is relevantly different since; what I'm searching for and still haven't found, is an actual explanation for what having a significant regression and no significant predictors means. How can this finding, if accurate, be explained? @Glen_b, I have seen that part of Jeromy's post (see my 2nd comment above) but how would I explain that in prose? Would it be accurate to say that the variables, together, have an effect on the outcome - but that their effects cannot be isolated from each other?


title: Not-significant F but a significant coefficient in multiple linear regression

text: I have a regression with two continuous predictors and one dichotomous predictor in Model 1 and two interactions of each of the continuous predictors with the dichotomous predictor in Model 2. The coefficient for one of the interaction terms is significant. However, the F test is not significant (neither for Model 1 nor for Model 2). Should I still interpret that significant interaction or should I just assume that neither the predictors, nor the interaction terms are useful in predicting the outcome variable?


title: Fitting multiple logistic functions at once with pymc

text: I have read a similar post here Fitting logistic function with pymc but it seems there is a rule that I shouldn't ask question in someone else's post. My approach is to fit 10 logistic functions to 10 different simulated data sets independently. The code written is done simultaneously with pymc's "assign_step_methods" function and sampling method is "pymc.AdaptiveMetropolis". There is one strange observations that needs your help. Every time when I run the script, there will be one data that doesn't fit the data at all. The estimated parameter seem like an outlier compared to the rest of the 9 set of estimated parameters that fit so well. The code seems to be bug free and able to fit 9 parameters. The priors are set within the range of the true parameters. The number of samples is set very large, almost 500,000. The only problem that I can think of is the likelihood function in which the "tau" used to in the Likelihood (Normal) function may be a little tight and it's a fixed value. I wonder if there are ways to check what went wrong with the fitting process. For example, there are lots of questions but how do I diagnose them ? a) If one of the 10 parameters explores in another space that results in rejections most of the time ? Is this related to the likelihood function ? Am I not understand the "Adaptive Metropolis" concept well ?   b) Why the algorithm keep missing one parameters ? I'm at the lost now. Your advice will be very helpful. Thank you. p.s. I'm currently using pymc2.2 and failed to install pymc3 to my Windows platform.


title: MCA vs Multinomial Logistic Regression

text: Lets say I have made survey using a sample of a given number of people, containing a set of 25 questions that have 6 possible answers (Fully agree/ Partially Agree/ Neutral/ Partially Disagree/ Fully disagree/ Don't know), and a few other questions concerning variables I'll be using to qualify them (such as gender, average income, region of residence, age, color, and so on). Now lets say I want to try and find correlations between gender/race/age/etc and the specific answers given to the set of the 25 questions. Lets say I also want to see if there might be correlations between answers. In other words, I want to find all possible correlations among these variables and the questions, and among the questions themselves. Should I be thinking about using a Multiple Correspondence Analysis for this objective, or a Multinomial Logistic Regression, or maybe both, or some other specific technique? I've been studying statistics by myself and I'm confronted with this problem and unsure of what path to follow. Both ways seem to want to find similar answers, but there seem to be nuances in choosing which one to use that I am missing.


title: ML for specific classification problem

text: I have a training dataset for classification problem $X \rightarrow y$. Where $X$ is an $n$th dimensional real vector, $y$ is an integer number in $\\{0, 1\\}$. I want to solve the next problem: predict as many lines where $y = 1$ as possible with error level no more than $\alpha = 10\%$. In other words, I want to split dataset ($N$ samples) in two sets: 1) $k$ - samples predicted to be $y = 1$, where $y = 0$ cases less than $k * \alpha$ 2) $N-k$ samples - I don't care about. I need to maximize $k$ with given $\alpha$ . I have a lot of experience in classification with ML, but this problem is rather unusual to me. Please help me how to formalize this problem for ML. Which and how should I use ML algorithm for the problem? I will be happy to get advice about general idea.


title: Why does fitdistr does not work with gamma?

text: I am tring to find a probability density for some waiting time, but I am having a hard time. `Fitdistr` does not work with Gamma. Am I missing something? Is there another density that fits better to this?               library(MASS)           wt = read.csv("http://isgeek.eu/SITN/waitingTime.csv", sep="\t", header = TRUE)     hist(wt$x, prob=TRUE)          # Estimate an exp proba     paraw <- fitdistr(wt$x,densfun="exponential")     curve(dexp(x, paraw$estimate[1]), 0,1300, add=TRUE, col="blue")     ks.test(wt$x, "pexp", paraw$estimate[1])           # Estimate a weilbull proba     paraw <- fitdistr(wt$x[wt$x!=0],densfun="weibull")     curve(dweibull(x, paraw$estimate[1], paraw$estimate[2]), 0,4000, add=TRUE, col="blue")     ks.test(wt$x, "pweibull", paraw$estimate[1], paraw$estimate[2])           # Estimate a gamma proba     paraw <- fitdistr(wt$x,densfun="gamma", list(shape = 1, rate = 0.1), lower = 40)     # this does not work...


title: Estimate power for linear model with Bernoulli-distributed error

text: I want to estimate a linear model for a phenotype $y_i$ defined by a liability threshold model, using observable data $x_i, c_i$. $$y_i = \mu + \beta x_i + \epsilon_i$$ $$\text{Pr}(\epsilon_i | x_i, \beta, \mu) \sim B(1, p_i) - p_i$$ $p_i$ is a complicated function of the normal CDF of the underlying error for the liability threshold model, which is itself based on a linear model with a different set of parameters, but for the purpose of this problem, I think the only important fact is that we cannot assume that $\epsilon_i$ is drawn from the same distribution for every individual $i$. I want to do maximum likelihood estimation of the parameters, and in particular, find the power of the estimates, but I'm not sure what to do because the errors are not only heteroskedastic, but depend on the individual.


title: When to scale/normalize for supervised learning algorithms?

text: I'm trying to understand which supervised learning algorithms require normalization/scaling of features. It appears that when an algorithm works by calculating the conditional probability (Naive Bayes, Linear Discriminant Analysis), then normalization is not necessary. However, for algorithms that use any kind of distance metric (K-nearest neighbors, Support Vector Machine), then normalization is necessary. Is this right? I'm a complete n00b so please jump in and correct me wherever appropriate.


title: What does it mean to build a so-called mean equation?

text: A series is denoted by Xt. How do we build a time-series model for Xt that is the "mean equation" of Xt? Is this another way of asking to build a stationary time-series?


title: What clustering algorithm can be used with a distance matrix and without feaures?

text: I have a dataset of binary files. I can't do feature extraction on them. I just computed the distance between every pair of file in the dataset with a distance metric (NCD = Normalized Compression Distance). So I have a distance matrix. My goal is to cluster these files. What is the best way to do that?


title: Show events on a trend plot

text: I am working on an algorithm to minimize a function by iterative improvement, and I am plotting the function's value for each iteration (pretty straightforward plot). I am feeding GnuPlot or PgfPlots an input file like:               "Solution 1"     2.22419e+007     2.22418e+007     2.22418e+007     ...      GnuPlot script:               set key autotitle columnheader     plot 'data.txt' with lines      PgfPlots script:               \begin{tikzpicture}     \begin{axis}[xlabel={Iteration}, ylabel={Solutions values}]          \addplot [blue] table[x expr=\coordindex, y index=0] {data.txt};          \end{axis}     \end{tikzpicture}      In my algorithm sometimes an event happens, namely I switch to another method for few iterations, or I have to relax some constraints because I can't get an improvement. I would like to show these events in the plot. I am looking for something like this (black vertical lines), but I am open to different alternatives as long as the meaning is the same: ![Ideal plot](http://i.stack.imgur.com/MSs3q.png) What is the best way to show events in this kind of plots? How to do it with GnuPlot and PgfPlots?


title: How to calculate the margin of error if the sample size is not given?

text: The mean and standard deviation of a random sample of $n$ measurements equal 33.9 and 3.3 respectively. Construct a 95% confidence interval for $n = 100$. Hence I am trying to calculate margin of error and then use that to derive the sample size and then to construct my diagram. But the information given seems not to be enough.


title: Why doesn't increasing N guarantee the sample mean will get closer to the population mean?

text: > When sampling from a population, the sample mean will always be closer to > the population mean as the sample size increases. Why is this statement wrong? My explanation is that as sample size increases, the standard error decreases and the sampling distribution will become normal. Am I right?


title: Question about classification with hidden Markov models using depmixS4

text: I am using the depmixS4 package to fit HMMs. I have three different classes of data and I have fitted 3 separate HMMs using the depmixS4 depmix and fit functions and given a new sequence of observations, I would like to be able to compute the probability (or log odds) of that new observation sequence given each of the three models that have already been trained. I understand this would normally be done using the Viterbi algorithm, but I have no idea how this can be implemented in depmixS4. Due to the functionality of depmixS4 I cannot use the HMM or HiddenMarkov packages Your help is greatly appreciated! Thank you in advance


title: Size of a test and level of significance

text: What is the difference between the two and why must the level of significance be always higher than or equal to the size of the test?


title: Some examples of errors in ESL?

text: I am reading the book _The Elements of Statistical Learning_ (Hastie _et al_ , chapter 7) and I am confused about the different kinds of errors mentioned in the book   1. Test error $$Err_\tau=E(L(y, \hat{f}(X)|\tau))$$   2. Expected prediction error $$Err=E(L(y, \hat{f}(X)))=E[Err_\tau]$$   3. Training error: $$\bar{err}=\frac{1}{N}\sum_{i=1}^{N}L(y_i,\hat{f}(x_i))$$   4. In sample error: $$Err_{in}=\frac{1}{N}\sum_{i=1}^{N}E_{Y^0}[L(Y_i^0,\hat{f}(x_i))]$$   5. Extra sample error: $$Err_\tau=E_{X^0,Y^0}(L(Y^0, \hat{f}(X^0)|\tau))$$ Can anyone give some idea to well understand these and some simple examples would be appreciated?


My binary logistic regression analysis (method=ENTER) model has all beta values are statistically non-significant. Should I include these predictors in the final binary logistic regression equation?


title: Test for whether a change is different in direction for subgroups?

text: I have a (large) population on which I randomly split into groups, and subjected one group to an experimental condition (call it change _A_ ) which demonstrated a significant improvement in a target variable overall. It happens that the type of the target variable is a proportion. However, the impact of change _A_ probably isn't the same for all subgroups. Suppose I split my population into men and women. Is change _A_ only better for women? I only care if the difference is a change in direction, that is, I don't care if it's 15% improvement for women and a 45% improvement for men, I only care if it's 10% worse for women and 20% better for men, in which case I should only apply the change for men. The change is costless but adds a little complexity to apply it only conditionally so if it's better for one subgroup and shows no significant difference for another subgroup, I'd rather apply the change across the board. I am only interested in whether the change is significantly positive for one group and significantly negative in another. The division into subgroups is somewhat arbitrary, men and women are presumably something on the order of a 50/50 split, but some of the segmentations I could make would be orders of magnitude smaller than the original population, like product brand usage or county or something like that, where the difference might not that significant for individual segments due the size of the subpopulation. Would it be a good idea to keep track of effectiveness of change _A_ , for individual counties, and only enable it for counties where the change was positive? In some counties the change was better, in some it was worse. How can I tell whether the differences are not due to random variation?


title: How to create mean score rounded to whole numbers in SPSS?

text: I have three variables measuring different aspects of performance each with a scale of 8 and I want to create a new Performance variable by computing these together. I have done this the following way: I've gone to "compute variable" in SPSS and chosen (var1+var2+var3)/3. Missing data is not a problem so I didn't have to add anything or that in the formula. Everything seems to be fine until this point and I only get full numbers between 1 and 8 by choosing "no decimals". However, when I take the descriptive statistics for this new variable, it gives me the same number multiple times. Because of this, I am not sure if I have created the new variable in the right way and if I can start running regression analyses


title: Measure of probability of obtaining a certain $\chi^2$ given the number of trials

text: I am calculating an equation with several random components a large number of times, and in order to select my final answer I choose that with the minimum value for $\chi^2$. I am hoping to find some way to decide how probable it is that I should end up with this value of $\chi^2_{min}$, to take into account the number of trials used. Initially I thought this was correct: $Probability = \frac{e^{-\chi^2_{min}}}{\sum_i^{N_{trials}}e^{-\chi^2_{i}}}$ So essentially the probability is the exponential of the lowest calculated $\chi^2$, divided by the sum of all calculated exp($-\chi^2$) values. This seems reasonable in some sense, because if there are a range of results with a low value for $\chi^2$ (as would be the case for a large number of trials) then the numerator would tend toward 1 whilst the demonimantor would grow with the trial number, making the overall probability that I select the best answer smaller. It doesn't really do what I would like it to do, however: in the exteme example where I have only 1 trial then this would always evaluate to '1', irrespective of the quality of the answer. In my equation I have several Gaussian random variables, and what I would really like is to have some measure of how 'far' a result is in terms of the means of these Gaussians, and include this in how I define the 'probability'. ie: if I have two results with similar $\chi^2$ but the random components of one are all very close to the means and the other a multiple $\sigma$ away, then I would select the first as the 'more probable' result. Is there any technique that exists to do something like this? Hopefully this is clear enough to comprehend - I'm pretty unsure of how next to proceed and any suggestions would be very welcome. **Editing to add a new idea:** I'm thinking of what I could do to take all of these elements into account is to do something like: $Probability = \frac{e^{-\chi^2_{min}} \times e^{-a} \times e^{-b}}{\sum_i^{N_{trials}}e^{-\chi^2_{i}} \times e^{-a} \times e^{-b}}$ where $a$ and $b$ (in $e^{-a}$ and $e^{-b}$) represent the distance from the mean of the Gaussian variables. Thoughts very welcome!


title: Is multiple logistic regression the right choice or should I use univariate logistic regression?

text: I have a set of data (~ 90 cases) and an outcome of a diagnostic test. I have collected factors that were determined before the test that could predict the outcome of the test. Now some of the data are binary, some are continuous (lab tests), one is categorical (the original diagnosis leading to the symptom with 5 categories). The statistician did (in SPSS) a multiple logistic regression on most of the parameters and then a backwards stepwise selection, which produced a model, one of the factors was strongly significant (<0.001), the others insignificant (two of them in the 0.01-0.05 range). The question from my colleagues was natural - why was this or that possible factor not included, and "these two factors are actually significant, or can we say nearing significance". I must say I don't like to interpret "nearing significance", but I have to work with my colleagues. I asked the statistician to perform another analysis, this time with the full set of factors, and the result is different, only one of the factors (the most significant one) stayed in the new model, is still strongly significant but there are other factors, also in the 0.01-0.05 range, which will bring me to the "nearing significance" problem again. I also tried with R and Rcmdr with the same data, and stepwise selection based on AIC or BIC produces different results (only the one factor remains the same). Now I see that given that various stepwise selection methods produce different models with insignificant factors, the presence of these factors in the model is just random. That's why I would like not to include them. Most of all, I don't want to interpret their presence in the model because the colleagues think that I should as they are "nearing significance". Question 1: To make things simpler, is it possible to run a "battery" of univariate logistic regressions on all the variables? This would solve both my problems - I would have one simple significant model with the most significant variable and wouldn't have to deal with "nearing significance" interpretation in the others. Question 2: Is there any correction of the significance to be used in multiple logistic regression as it is with multiple comparisons, so I can deal with the "near significance" argumentation?


title: Bayesian Perceptron - how can I generate many different perceptrons?

text: I am going to implement the Bayesian version of a perceptron that I read in the Statistical Mechanics of learning, by Engel-Van Den Broeck. The idea to improve the performance is to use many Gibbs perceptrons J to define a golden one using their center of mass $J^\beta$. Formally: $$J^\beta = \lim_{m\rightarrow \infty} \sqrt N \frac{\sum^M_{i=1}J^{(i)}}{\sqrt{(\sum^M_{i=1}J^{(i)})^2}}.$$ The problem is: how should I compute these many $J$? Should I just use many different initial condition? Should I use different $J$ taken during different epochs of the perceptron learning (this in particular is my thought)? Should I perturb data? Should I perturb one specific $J$ in many different random ways?


title: GB-RBM unable to learn and generate simple 2D-Motions?

text: I am trying to apply a GB-RBM to a variation of the _LASA Handwriting Dataset_ to be able to generate new examples. My dataset contains motion trajectories of simple shapes like 'S' shapes, spirals, 'C' shapes etc. Each motion is represented by a list of $n$ 2D-points $(x_i, y_i)_n$. The dataset is normalized, so that all trajectories have the same starting and end point which are $(-1,0)$ and $(0,0)$. This is what a trajectory could look like: ![Spiral motion from the LASA Handwriting Dataset](http://i.stack.imgur.com/rMjDw.png) I merge each of these trajectories into one vector by simply appending all $x$ and $y$ components: $v_k = (x_1 \dots x_n, y_1 \dots y_n)$. I trained a GB-RBM on this data with many different learning rates, different hidden layer sizes, with learned and fixed variance... Somehow the RBM is unable to generate all kinds of shapes that occur in the dataset. The generation process almost always generates one distinct shape (spiral) - already after a handful of steps of sampling (starting from a random-vector). Are the gaussian linear units fundamentally unable to learn this kind of distribution properly? Or is this a problem with the binary hidden units? Or something completely different?


title: Conditional or Joint Probability under Various distributions

text: In various statistical models the baseline equation (like in Naive Bayes $$\mathrm{classify}(f_1,\dots,f_n) = \underset{c}{\operatorname{argmax}} \ p(C=c) \displaystyle\prod_{i=1}^n p(F_i=f_i\vert C=c)$$ changes under various distributions (like for Gaussian Naive Bayes $$p(x=v|c)=\frac{1}{\sqrt{2\pi\sigma^2_c}}\,e^{ -\frac{(v-\mu_c)^2}{2\sigma^2_c} } ,$$ again for Multinomial Gaussian Naive Bayes, it is, $$p(F|C) = \frac{(\sum_i F_i)!}{\prod_i F_i !} \prod_i {p_i}^{F_i} $$ Similar constructions are there for Hidden Markov Model, and I believe each statistical model may have same issue. I tried to check my Statistics book over the definitions of respective distributions, but the answer is bit presumed and may not be satisfactory, so I thought to ask the experts over here. I am slightly confused how these equations are changing like what is the parameter of change? Is it if I take the base equations and apply appropriate distribution models, I would get the respective equations? If any one of the esteemed members may please help.


title: Combining data from different sources

text: I want to combine data from different sources. Let's say I want to estimate a chemical property (e.g. a partitioning coefficient): I have some empirical data, varying due to measurement error around the mean. And, secondly, I have a model predicting an estimate from other information (the model has also some uncertainty). How can I combine those two datasets? [The combined estimate will be used in another model as predictor]. Meta-analysis and bayesian methods seem to be suitable. However, haven't found many references and ideas how to implement it (I am using R, but also familiar with python and C++). Thanks. **Update** Ok, here's a more real example: To estimate the toxicity of a chemical (typically expressed as $LC_{50}$ = concentration where 50% of animals die) lab-experiments are conducted. Happily the results of the experiments are gathered in a database (EPA). Here are some values for the insecticide Lindane:               ### Toxicity of Lindane in ug/L     epa <- c(850 ,6300 ,6500 ,8000, 1990 ,516, 6442 ,1870, 1870, 2000 ,250 ,62000,              2600,1000,485,1190,1790,390,1790,750000,1000,800     )     hist(log10(epa))          # or in mol / L     # molecular weight of Lindane     mw = 290.83 # [g/mol]     hist(log10(epa/ (mw * 1000000)))      However, there are also some models available to predict toxicity from chemical properties (QSAR). One of these models predicts toxicity from the octanol/water partition coefficient ($log~K_{OW}$): $$ log~LC_{50} [mol/L] = 0.94~(\pm 0.03)~log~K_{OW}~-~1.33 (\pm~0.1)$$ The partitioning coefficient of Lindane is $log~K_{OW} = 3.8$ and the predicted toxicity is $ log~LC_{50} [mol/L] = -4.902 $.               lkow = 3.8     mod1 <- -0.94 * lkow - 1.33     mod1      Is there a nice way to combine these two different informations (lab experimens and model predictions)?               hist(log10(epa/ (mw * 1000000)))     abline(v = mod1, col = 'steelblue')      The combined $LC_{50}$ will be used later on in a model as predictor. Therefore, a single (combined) value would be a simple solution. However, a distribution might be also handy - if this is possible in modelling (how?).


title: Tractability of mutual information-augmented ensemble classification algorithms

text: I am seeking to augment random forest classification using Shannon-Weaver mutual information as a metaheuristic to partition candidate datasets. Specifically, I am trying to determine if such an approach is both tractable and offers an improved convergence time over bruteforce methods and random forest classification on its own without a major reduction in classification. Suppose we have a collection of data sets, $K$, that contains rows $\\{k1, k2, ..., kn\\}$. Each of these rows carries an unknown language $L$ in the context of an unknown, uniformly-applied noise function. Our hypothesis is that at least some subset of these data rows share information: that is, for some rows, $I(kx,ky) > 0$ for ${kx,ky} \in K$. Additionally, we hypothesize that multivariate mutual information $I(kx,ky,...) \not= 0$ for some set of multivariate relations. That is, that our graph contains compositions that are more informative in the context of multiply-chained relationships, or contain some XOR-shaped relationship. If we attempt to bruteforce this space, because mutual information is commutative, we require n-binomial-n time and space. That is to say, to compute the mutual information of all combinations within collection $D$, we require $\sum\limits_{x=1}^n\binom{n}{x}$ time, assuming for simplicity that the comparison of each row is an $O(1)$ operation. We would then prune output data set using a fixed cutoff value, $v$, retaining only connections where $I(kx,...) > v$. We would then sort this set in descending order of the number of terms and use this to construct an undirected graph, by connecting features that share positive mutual information and forming branches when multivariate mutual information is negative.    The obvious reduction to this complexity, at an error-bound cost in robustness, is to use random forest analysis. I wish to simplify this computation further, by computing the pairspace mutual information for each pair of rows, then using a cutoff value similar to $v$, construct sets that only contain the participating nodes of each of these graphs. These sets would then be fed into our random forest classifier to determine the shape of the graphs, at an expected reduced time for convergence. Is this approach tractable? Even if so, are there better approaches to the pruning stage of this problem? (Originally posted on CSTheory.SE)


title: Proper interpretation for error analysis

text: Given an arbitrary discrete distribution and an observed distribution coming from a Monte Carlo simulation, my goal is to be able to say whether or not the observed distribution is the same is the given distribution such that I correctly identify a distribution which is different with probability greater than 1 in 10,000. I should say that the Monte Carlo simulation is run at 10^10 iterations. So far, I have been using a combination of confidence intervals (derived from the given distribution) to test the observed mean as well as the chi-square test. Originally, it seemed to me that this combination would provide, over multiple iterations, the desired precision for which I am aiming. Upon further thought, however, it has occurred to me that these tests are not independent since I am using them on the same observation. Consequently, I can only do as well as the most precise test. Is this true? I have several other questions: Is my interpretation of confidence intervals correct when I say, when using a 95% confidence interval (for instance), that the probability that the observed mean falls outside of the confidence interval when the observed distribution is the same as the given distribution is 5%? Similarly, if the above is correct, is the interpretation of the chi-square result similar? In other words, with an alpha value of 0.05, is it true that the probability that the chi-square result falls outside of this bound when the observed distribution is the same as the given distribution is 10%? Finally, is there a good method for achieving my desired precision? Thank you in advance!


title: How does the number of parameters scale with the number of features?

text: I am having trouble wrapping my brains around this seemingly simple, yet conceptual problem. As the title suggests, the question is: How does the number of parameters required for Naïve Bayes, LDA, and Logistic Regression (unregularized) scale as a function of the number of features? If Y is boolean with $X_1 ... X_p$ features: I understand that for logistic regression we have a linear relationship between the number of parameters as based on the number of features based on the degrees of freedom in a null model subtracted from the residual degrees of freedom when a glm function is performed. i.e. $1 + p$ For Naive Bayes, I am think I understand that, but I am not certain: $P(Y) = p \cdot (\mu(Y=1), \mu(Y=0), \sigma(Y=1), \sigma(Y=0)) = 1 + 4p$ For LDA? I'm lost. Does anyone have any insight?


title: Do you need causal models when doing counterfactual predictions?

text: I am modeling the impact the number of a certain type of company (bottom of pyramid (BOP) companies, ie. companies that cater to the poorest consumers) have on market price. I considered the following two models: (1) The first model was designed to capture causal effects; we instrumented the number of BOP companies to estimate the causal relationship between the number of BOP companies on average market price. I happy with this model. (2) I am a bit puzzled about the second model I am considering. What I am trying to do now is to estimate the total income had no BOP company entered the market, and the total market income had these companies entered the market earlier. There are around 100 different markets being considering, covering around 10 years. Suppose we were to estimate the total income of these markets if no BOP company entered in any of them. Using the model from part (1), and inputing the number of BOP companies as zero, I can predict what would have happend if there were no BOP companies over the observed years. From this point, I can basically compute the total market revenue. My question is: QUESTION: Is it right to use the model of step (1) on step (2)? Do I really have to care about endogeneity when what I want is to do a prediction of the type in (2)? More troublesome, it is actually wrong to restrict myself to a model of type (1) when I want to do step (2)?


title: Single-class vs two-classes hypothesis testing

text: Consider a Gaussian mixture model of unknown parameters, having just one or two components. I would like to design a statistical test to decide the number of classes ($H_0$: single component model vs. $H_1$: two components model). I don't know how to formulate this. Can someone help me get started?


title: Impose several restrictions on systems of eqns in R using systemfit package

text: I'm estimating systems of simultaneous eqns in R, and trying to impose homogeneity and symmetry restrictions on the equations. I'm using the "systemfit" package. I've managed to impose one restriction on one of the eqns using 'restrict.matrix' (as you see from the example below). How do I impose all restrictions at the same time (restrict1, restrict2, restrict3 and so on)?               eqn1 <- s1 ~ p1 + p2 + p3 + p4 + q     eqn2 <- s2 ~ p1 + p2 + p3 + p4 + q     eqn3 <- s3 ~ p1 + p2 + p3 + p4 + q      eqlist <- list (eqn1, eqn2, eqn3)          system_eqn <- systemfit(eqlist, method="SUR")          restrict1 <- "eq1_p1 + eq1_p1 + eq1_p3 + eq1_p4 = 0"     restrict2 <- "eq2_p1 + eq2_p1 + eq2_p3 + eq2_p4 = 0"     restrict3 <- "eq1_p2 - eq2_p1 = 0"          system_eqn_restrict <-(eqlist, method="SUR", restrict.matrix = restrict1)


title: Whether the population exists in an experimental study?

text: In my oral exam one of my teachers asked whether the population exists in an `experimental study`. I answered: "Yes. the population exists in an `experimental study`. And I will get it from my interested subject." But he rejected the answer and told me that the population does not exist in an `experimental study` but exists in an `observational study`. I have not understood it. If the population does not exist in an `experimental study`, then how can I get my sample ?


title: What is a second order stationary process?

text: I was wondering how his "second-order stationary process" is defined in Brockwell and Davis' Introduction to Time Series and Forecasting: > The class of linear time series models, which includes the class of > autoregressive moving-average (ARMA) models, provides a general framework > for studying stationary processes. In fact, every **second-order stationary > process** is either a linear process or can be transformed to a linear > process by subtracting a deterministic com- ponent. This result is known as > Wold's decomposition and is discussed in Section 2.6. In Wikipedia, > The case of second-order stationarity arises when the requirements of strict > stationarity are only applied to pairs of random variables from the time- > series. But I think the book has a different definition from Wikipedia's, because the book uses stationarity short for wide-sense stationarity, while Wikipedia uses stationarity short for strict stationarity. Thanks and regards!


title: How to access coefficient from VGAM package (generalized Poisson regression)

text: I am interested to know, how to get the coefficient from VGAM package (generalized Poisson regression), because, when I want to get the coefficient, is just can to get the estimate with @ to get the estimate. An other coefficient can to get all of the part from coefficient with "$"... My program like this:               library("VGAM")       fit <- vglm(y ~ x1, genpoisson,trace = TRUE)       summary(fit)       hasilgp <- summary(fit)       koefgp <- hasilgp@coefficients        And only the estimate can to attach... Any one help me...


title: Autocorrelation and evidence of iid

text: Suppose I have the first seven autocorrelations for some variable $x$.   And suppose they are -0.2, 0.15, -0.05, -0.10, -0.05, -0.14, 0.04 How can this be used as evidence of my data being or not being IID (suppose $n=150)$.


title: Distribution of sum of squares error for linear regression?

text: I know that distribution of sample variance   $$ \sum\frac{(X_i-\bar{X})^2}{\sigma^2}\sim \chi^2_{(n-1)} $$ $$ \sum\frac{(X_i-\bar{X})^2}{n-1}\sim \frac{\sigma^2}{n-1}\chi^2_{(n-1)} $$ It's from the fact that $(X-\bar{X})^2$ can be expressed in matrix form, $xAx'$ (where A: symmetric), and it could be again be expressed in: $x'QDQ'x$ (where Q: orthonormal, D:diagonal matrix). What about $\sum(Y_i-\hat{\beta}_0-\hat{\beta}_1X_i)^2$, given the assumption $(Y - \beta_0 - \beta_1X)\sim \mathcal{N}(0, \sigma^2)$ ? I figure $$\sum\frac{(Y_i-\hat{\beta}_0-\hat{\beta}_1X_i)^2}{\sigma^2}\sim \chi^2_{(n-2)}.$$ But I have no clue how to prove it or show it. Is it distributed exactly as $\chi^2_{(n-2)}$?


title: Is it appropriate to run correlations first and then a regression?

text: In my study I have several predictor variables as well as dependent variables/outcomes. I don't have much theoretical pressure to choose any variable over another, making decisions about regressions frustrating. So I ran the correlations between all variables, IVs and DVs to see which had relationships. Is this frowned upon? Can I use this information to inform where to look further? Are correlations susceptible to familywise error?


title: Problem with creating proper linear regression model for my experiment

text: I am working with microarray data and I am trying to find differentially expressed genes affected by gut microbiota in mice. I have 6 subjects, out of which I got 4 different samples one from L-negative cells and one from L-positive cells from both colonic and ileal tissues. In the following you can see the designed matrix for each microarray file: > >     file    subject   gutMic    tissue    Lcells >     -------------------------------------------- >     IPCR1       1     CR        ileum     Lpos >     CPCR1       1     CR        colon     Lpos >     INCR1       1     CR        ileum     Lneg >     CNCR1       1     CR        colon     Lneg >     IPCR2       2     CR        ileum     Lpos >     CPCR2       2     CR        colon     Lpos >     INCR2       2     CR        ileum     Lneg >     CNCR2       2     CR        colon     Lneg >     IPCR3       3     CR        ileum     Lpos >     CPCR3       3     CR        colon     Lpos >     INCR3       3     CR        ileum     Lneg >     CNCR3       3     CR        colon     Lneg >     IPGF1       4     GF        ileum     Lpos >     CPGF1       4     GF        colon     Lpos >     INGF1       4     GF        ileum     Lneg >     CNGF1       4     GF        colon     Lneg >     IPGF2       5     GF        ileum     Lpos >     CPGF2       5     GF        colon     Lpos >     INGF2       5     GF        ileum     Lneg >     CNGF2       5     GF        colon     Lneg >     IPGF3       6     GF        ileum     Lpos >     CPGF3       6     GF        colon     Lpos >     INGF3       6     GF        ileum     Lneg >     CNGF3       6     GF        colon     Lneg >   Now I would like to analyse the effect of gut microbiota in L-negative comparing to L-positive cells regardless of tissue type (i.e. GF ~ CR in Lneg and compare it with GF ~ CR in Lpos). I know for data coming from same subject I should consider correlation raising from measurement in the same subject to reduce biased result. But, here I don't know how I can reduce effect of samples coming from different tissues of one subject?!! As an example files INGF3 and CNGF3 are coming from same subject but different tissues and I would like to sum them up as consider them replicates of one sample.


title: Is is valid to calculate the mean of multiple neural networks results?

text: Recently I've been playing around with Artificial Neural Networks models and I came across a problem that I can't figure out. Essentially, I'm playing around with a model trained on a relatively small sample from which the forecast performance of the trained neural networks model is limited. However, I noticed that if I run the model multiple times, the average of the sample forecast tends to be more satisfying. My question is, is this a valid way to utilize the forecast results? Is there some other way around that I should seek to achieve a similar goal?


title: Is the survival function the same as the upper tail?

text: Are the two synonymous? The reason I am asking is that I have a paper that calculates a certain p-value as the upper tail of the hypergeometric distribution: $\Sigma_{k}^{m} = \frac{\binom{m}{k}\binom{N-m}{n-k}}{\binom{N}{n}}$ I want to implement their proposal using Python. Scipy offers a hypergeometric distribution and it can calculate the survival function of that distribution from a given point. So I am wondering, are the two the same? Can I just use that function?


title: Graphing low-end prices of goods in a market over time

text: I have a large sample of market data - that is the prices and amount of goods being sold by vendors at specific, but inconsistent, points in time. Assuming I am a buyer and the quality of the goods are constant across the entire market (consumable goods by one manufacturer, let's say), how would you graph this data? Keep in mind that outliers may exist (fire sales or inflated asking prices). Seeing as I'm only concerned with the lowest-priced goods at a specific moment in time, a market average wouldn't do much for me unless I was buying most of the market's supply. I considered a box plot, but that also tends to be mean/median-biased. A simple line/area/bar graph of the lowest prices over time doesn't seem to give enough information, especially if there isn't much volume available at the lowest price. *I'm not a stats guy, so go easy on me!


title: Bounded in probability and finite expectation

text: Let $x_t = O_p(1)$, meaning that for all $\varepsilon > 0$ there exists $M_{\varepsilon} < \infty$ s.t. $P(|X_t| > M_{\varepsilon}) < \epsilon$ for all $t \in \mathbb{N}$. Does it imply that $\sup_{t} Ex_t < \infty$, namely, whether the expectation of a random variable which is bounded in probability is finite. If it is not the case, could you give me a counter example?


title: R code for Kutner et al's Applied Linear Statistical Models?

text: I am self-studying Applied Linear Statistical Models by Michael Kutner, Christopher Nachtsheim, John Neter, William Li. Is there R code for the methods in the book?


title: Measuring goodness-of-fit in a model that combines two distributions

text: I have data with a double peak that I'm trying to model, and there's enough overlap between the peaks that I can't treat them independently. A histogram of the data might look something like this: ![alt text](http://i.stack.imgur.com/EuTFp.png) I've created two models for this: one uses two Poisson distributions, and the other uses two negative binomial distributions (to account for overdispersion). What's the appropriate way to tell which model fits the data more accurately? My initial thought is that I could use a Kolmogorov-Smirnov test to compare each model to the data, then do a likelihood ratio test to see if one is a significantly better fit. Does this make sense? If so, I'm not exactly sure how to perform the likelihood ratio test. Is chi-squared appropriate, and how many degrees of freedom do I have? If it helps, some (very simplified) R code for the models might look something like this:               ## inital data points     a <- read.table("data")          #create model data     model.pois = c(rpois(1000000,200),rpois(500000,250))     model.nb = c(rnbinom(1000000,200,0.5),rnbinom(500000,275,0.5)          #Kolmogorov-Smirnov test     #use ks.boot, since it's count data that may contain duplicate values     kpois = ks.boot(model.pois,a)     knb = ks.boot(model.nb,a)          #here's where I'd do some sort of likelihood ratio test     # . . .      * * * Edit: Here's an image that may explain the data and the distributions I'm fitting better. It's totally clear from the visualization that the second model (using the negative binomial dist to account for overdispersion) is a better fit. I'd like to show this quantitatively, though. ![alt text](http://i.stack.imgur.com/NE6Vq.png) (red - data, green - model)


title: Log likelihood improves with addition of a nonsignificant variable

text: In this question I asked about changes to AIC when adding a variable. It turns out to be partly due to the way SAS figures AIC. However, I now have two models where the log likelihood improves a lot: This code               title 'Test - only forced';     ods select FitStatistics SolutionF;     proc mixed data = menfat.data method = ml;      class &classvar;      model  bmd_legneck_change_per_year = &bmddef bmdlegneck_1/solution;     run;      had -2LL of 1631.9 adding one variable with this code               title 'Test - add packyrs';     ods select FitStatistics SolutionF;     proc mixed data = menfat.data method = ml;      class &classvar;      model  bmd_legneck_change_per_year = &bmddef packyrs bmdlegneck_1/solution;     run;      had a -2LL of 1607.6 (which is highly significant difference of 24.3 on 1 df) but the parameter estimate for packyrs was 0.0016 with an SE of 0.005 and a p of 0.77. Could this be due to improvement of the fit of other variables in the model? EDIT: Responding to Macro; the model is a linear regression where the dependent variable is a change per year of a physical quality (bone density) and the independent variables are a mix of yes/no variables and continuous ones.


title: Will the silhouette formula change depending on the distance metric?

text: I am using Silhouette width to compute the best value for k in k-means. As I am performing document clustering, I am calculating the values of `a` and `b` as follows in Python:               a = distance(data[index], centroids[clusters[index]], metric=metric, p=p)     b = min([distance(data[index], c) for i,c in enumerate(centroids) if i != currentindex])     score = float(b - a) / max(a, b) if max(a, b) > 0 else 0.0      The original formula from the Wikipedia page is the following: ![enter image description here](http://i.stack.imgur.com/2mdKh.png) I am using a `cosine` similarity measure to compute the distances and was wondering if this formula needs to be changed or it can be left as is. In the above snippet of code, the function `distance` computes the cosine similarity metric.


title: Test for independence of random variables

text: I have a time series of data (about 300-750 elements, depending on the sample) and a model that has some random residues. I used the Kolmogorov–Smirnov test to make sure that the normality hypothesis can't be rejected, so I assume that the residuals are normally distributed. But now I guess I should test if they are independent of each other - so that there is no autoregression? Which test should I use (preferably one that is easily implementable in java)?


title: alpha-quantile of chi square

text: Please take a look first at this link: http://en.wikipedia.org/wiki/Ljung%E2%80%93Box_test#Formal_definition It is written, $\chi_{1-\alpha,h}^2$ is the $\alpha$-quantile of the chi- squared distribution with $h$ degrees of freedom. Does $\chi_{1-\alpha,h}^2$ mean:   1. The area below the graph of $\chi^2(h)$ from $\chi_{1-\alpha,h}^2$ to $\infty$ is $\alpha$, or   2. The area below the graph of $\chi^2(h)$ from $\chi_{1-\alpha,h}^2$ to $\infty$ is $1-\alpha$?


title: Check if a group of signals have a significantly lower or higher value from another group using GLM

text: I have a question that is a kinda puzzling me and hope that I can find some answers here :). Suppose I have many signals, let's say four $X_1, X_2, X_3, X_4$ Suppose I have a set of windows $w_1, w_2, \dots , w_n $, in each window each signal follow a Poisson distribution, i.e: $ X_i \sim Poisson(\lambda_i) $. What I want to do is to find windows that show a high value for some signals and not for others, for example they have high value for $X_1, X_2$ and low signal for $X_3, X_4$ One idea that come to my mind is to use an ANOVA-like method in which the first group contains $ X_1, X_2$ and the second contains $ X_3, X_4$ and test if the values of group are significantly **bigger** or **lower** than the other. Because the data is not normally distributed I though of using a GLM model how ever it seems that existing glm package in R and python test for equality (hope to be wrong). Does this approach seem reasonable ? if so how to perform this test using R or python statsmodel package. Thanks in advance,


How can I improve the accuracy of my logistic regression code, which tests the accuracy using the 10-fold cross-validation technique? I have implemented this code using `glmfit` and `glmval`. The desired accuracy is somewhat higher and it requires the parameters to be found using maximum likelihood estimator. Also, when I run this code in MATLAB, I get the following error > Warning: X is ill conditioned, or the model is overparameterized, and some > coefficients are not identifiable. You should use caution in making > predictions. In glmfit at 245 In LR at 8 The code is:               function LR( X,y)     y(y==-1)=0;     X=[ones(size(X,1),1) X];     disp(size(X,2));     indices = crossvalind('Kfold',y,10);     for i = 1:10         test = (indices == i); train = ~test;         b = glmfit(X(train,:),y(train),'binomial','logit');         y_hat= glmval(b,X(test,:),'logit');         y_true=y(test,:);         error(i)=mean(abs(y_true-y_hat));     end     accuracy=(1-error)*100;     fprintf('accuracy= %f +- %f\n',mean(accuracy),std(accuracy));     end


title: Writing an ordinal (proportion odds model) as function of binary regression

text: I wonder if it is possible to write a proportional odds model as some function of binary logistic models. Indeed I have data on an ordinal response scale (mild, moderate, poor). There are options to dichotomous the data because, we will only need to apply two treatments. But of-course an ordinal model can be use for the classification, then the choice of the two treatments can be applied accordingly. The question is which model does well on the data? How do we quantify "do well"? Available solutions are: 1) Build the models and look at the AUC for binary and VUS for ordinal and compare them. 2)Look at the Brier score for ordinal and binary model. But I still think the fundamental solutions lie on the effects we are modelling ($\beta$) how do they relate for the two models.


title: Is this identifiable?

text: I am interested in the following model : ($1 \leq i \leq p$, $1 \leq j \leq n_{i}$) $$y_{i,j} = A (1+a_{i})(t_{i,j}+\gamma_{i}) + \varepsilon_{i,j}$$ where $A \in \mathbb{R}$, $(a_{i})_{1 \leq i \leq p}$ is a $p$-random sample from $\mathcal{N}(0,\sigma_{a}^{2})$ and $(\gamma_{i})_{1 \leq i \leq p}$ is a $p$-random sample from the distribution $\mathcal{N}(0,\sigma_{\gamma}^{2})$, independent from $(a_{i})_{1 \leq i \leq p}$. $(\varepsilon_{i,j})_{i,j}$ is a random sample from $\mathcal{N}(0,\sigma^{2})$, independent from $(a_{i})_{i}$ and $(\gamma_{i})_{i}$. Given some observations $(y_{i,j},t_{i,j})_{i,j}$, I would like to estimate the model parameters ($A,\sigma_{a}^{2},\sigma_{\gamma}^{2},\sigma^{2})$. I am wondering whether this is possible. I created some synthetical data using the model in Matlab and tried to estimate the parameters of the model using `nlmefit`(for nonlinear mixed-effects models) but the estimated values for $A$ and $\sigma_{a}$ were (respectively) three times smaller (and three times) bigger than the ones I used to create the data. This makes me wonder whether the model is identifiable. Can you see another reason to this ?


title: KL divergence between an uninformative (?) Gaussian and a Gaussian

text: I have to calculate the KL divergence between a distribution $q$ and a prior distribution $p$, both of which are univariate Gaussians, i.e. $KL(q|p), q \sim \mathcal{N}(\mu, \sigma^2), p \sim \mathcal{N}(\mu', \sigma'^2)$. This term is part of a larger formula, which is justified in some way not relevant to this question. Now, to be honest, I don't want to put _any_ asusmptions on $p$ except that it is Gaussian. My intuition is to just say that if I do that, I can just say $p=q$ and thus $KL(p|q) = 0$. I wonder if there is a way to phrase this intuition into proper math. I read about non-informative priors, but found nothing about calculating KL divergences in this scenario. **Update:** I try to make the question more clear. I have two distributions. One is the prior, the other driven by data. My prior is not tied to any parameter ranges (e.g. in form of a conjugate prior), but only specified in its functional form (i.e. disitrbution family). In typical Bayesian frameworks, such a thing is called an informative prior. Does the same concept exist for KL based objective functions?


title: How can I continue training a neural net in R?

text: I have data set that is up-dated every month. I've already trained a neural network on a development sample and validated it. For the building process I used R's nnet package. I also tried neuralnet. I am happy with the model. However every month I will get an up-dated dataset. How can I continue to train the existing network using either of those packages? Thank you very much in advance, N


title: Hyper-geometric Probability Distribution in item defect problem

text: I was looking at the following example but couldn't get any solution... > Items from a large lot are examined one by one until r items with a rare > manufacturing defect are found . The proportion of items with this type of > defect in the lot is known to be $'p'$. Let $X$ denote the number of items > needed to be examined. Derive the probability distribution of $X$, and find > $E\left (X \right )$. I think this is related to hyper-geometric function but couldn't go further.


title: Why is the conditional expectation of prediction error in regression not zero?

text: The conditional expectation of the error in regression is: $E[Y-X\beta|X=x]$ is not equal to 0. Why is this the case? If you fix all the predictor variables, why does $E[Y]$ - $X\beta$ not equal to 0?


title: How to compare different distributions with reference truth value in Matlab?

text: I have production (q) values from 4 different methods stored in the 4 matrices. Each of the 4 matrices contains q values from a different method as:               Matrix_1 = 1 row x 20 column           Matrix_2 = 100 rows x 20 columns           Matrix_3 = 100 rows x 20 columns           Matrix_4 = 100 rows x 20 columns       The number of columns indicate the number of years. 1 row would contain the production values corresponding to the 20 years. Other 99 rows for matrix 2, 3 and 4 are just the different realizations (or simulation runs). So basically the other 99 rows for matrix 2,3 and 4 are repeat cases (but not with exact values because of random numbers). Consider `Matrix_1` as the reference truth (or base case ). Now I want to compare the other 3 matrices with `Matrix_1` to see which one among those three matrices (each with 100 repeats) compares best, or closely imitates, with `Matrix_1`. How can this be done in Matlab? I know, manually, that we use confidence interval (CI) by plotting the `mean of Matrix_1`, and drawing each distribution of `mean of Matrix_2`, `mean of Matrix_3` and `mean of Matrix_4`. The largest CI among matrix 2, 3 and 4 which contains the reference truth (or `mean of Matrix_1`) will be the answer.               mean of Matrix_1 = (1 row x 1 column)          mean of Matrix_2 = (100 rows x 1 column)          mean of Matrix_3 = (100 rows x 1 column)          mean of Matrix_4 = (100 rows x 1 column)      I hope the question is clear and relevant to SA. Otherwise please feel free to edit/suggest anything in question. Thanks!


title: What's an intuitive way to explain SS decomposition for a Mixed ANOVA with 1 B/S and 1 W/S variable?

text: For teaching I need to explain the following diagram to undergraduate students. The diagram is taken from the 7th Edition of Howell's "Statistical Methods for Psychology". The students are familiar with how to calculate the Total SS, Between SS, and Within SS for a One-Way Independent groups design, but don't know anything more than that. ![enter image description here](http://i.stack.imgur.com/EL0eM.png)


title: Are "random sample" and "iid random variable" synonyms?

text: I have been facing hard time understanding meaning of "random sample" as well as "iid random variable". I tried to find out the meaning from several sources, but just got more and more confused. I am posting here what I tried and got to know: Degroot's Probability & Statistics says: > Random Samples / i.i.d. / Sample Size : Consider a given probability > distribution on the real line that can be represented by either a p.f. or a > p.d.f. $f$. It is said that $n$ random variables $X_1 , . . . , X_n$ form a > random sample from this distribution if these random variables are > independent and the marginal p.f. or p.d.f. of each of them is $f$. Such > random variables are also said to be independent and identically > distributed, abbreviated i.i.d. We refer to the number n of random variables > as the sample size. But one of the other statistics book I have says: > In a Random Sampling, we guarantee that every individual unit in the > population gets an equal chance(probability) of being selected. So, I have a feeling that i.i.d.s are elements that construct random sample, and the procedure to have random sample is random sampling. Am I right? P.S.: I am very confused about this topic, so I will appreciate elaborate reply. Thanks.


title: Conjugate prior for a Gamma distribution

text: I need to update the failure rate (given as deterministic) based on new rate of failure about the same system (it is a deterministic one too). I read about conjugate priors and Gamma distribution as a conjugate for the Poisson process. Also, I can equate the mean value of Gamma dist. ($\beta/\alpha$) to the new rate (as a mean value) but I do not have any other information such as standard deviation, Coefficient of Variation, 90th percentile value,...etc. Is there a magic way to manipulate that and find parameters for the prior Gamma hence I get the posterior which Gamma too?


title: Can someone explain the C-Index in the context of hierarchical clustering?

text: This is a followup to this question. I am currently trying to implement the C-Index in order to find a near-optimal number of clusters from a hierarchy of clusters. I do this by calculating the C-Index for every step of the (agglomerative) hierarchical clustering. The problem is that the C-Index is minimal (0 to be exact) for very degenerated clusterings. Consider this: $c = \frac{S-S_{min}}{S_{max}-S_{min}}$ In this case $S$ is the sum of all distances between pairs of observations in the same cluster over all clusters. Let $n$ be the number of these pairs. $S_{min}$ and $S_{max}$ are the sums of $n$ lowest/highest distances across all pairs of observations. In the first step of the hierarchical clustering, the two closest observations (minimal distance) are merged into a cluster. Let $d$ be the distance between these observations. Now there is one pair of observations in the same cluster, so $n=1$ (all other clusters are singletons). Consequently $S=d$. The problem is that $S_{min}$ also equals $d$, because $d$ is the smallest distance (that is why the observations where merged first). So for this case, the C-Index is always 0. It stays 0 as long as only singleton clusters are merged. This means the optimal clustering according the C-Index would always consist of a bunch of clusters containing two observations, and the rest singletons. Does this mean that the C-Index is not applicable to hierarchical clustering? Am I doing something wrong? I have searched a lot, but could not find any suitable explanation. Can someone refer me to some resource that is freely available on the internet? Or, if not, at least a book I may try to get at my universities library? Thanks in advance!


title: Can confidence intervals be added?

text: For example, say we have three coins. Due to several tests, we have 95% confidence intervals for X (the heads percentage) of each coin. Coin A is 0-5%, coin B is 1-9%, and coin C is 40-70%. [For the curious, these CI's were created using inverse binomial CDF, based a small sample of ~30 flips of the coin. Not surprisingly, we have no clue about the composition of the coin-bag these coins were drawn from.] Can we say anything about the number of heads we would expect to see if we flipped all three? Or if we chose one coin randomly and flipped it? Based on the computed CI's alone? What if we can also use the info that led us to those CI's? **Edit: let me tell you much more about my problem.** We have 50 coins, all of which we suspect are very low probability heads (well under 5%). We do 30 trials, where each trial consists of flipping all the coins and recording all the results. The value we seek is a CI on the average number of heads in a trial. One approach would be to record the 30 "numbers of heads", and do a standard normal approximation. The problem is that many of these 30 numbers are 0 and 1, so while we get the correct mean, I fear the resulting CI is way tighter than it has any right to be. What I would like to do is calculate a CI for each coin, using it's 30 binary outcomes. I have calculated that for 0 heads out of 30, this leads to a 95% CI of 0-11.5% (with a mean of 0, of course). I did this by substituting values for P in "cumulative distribution function(1) of binomial distribution n=30 p=P" in Wolfram Alpha until the result was around 0.025. For 1 heads out of 30, I get a CI of 0-17.2% Once I have done this, can I add up the CI's in any way?


title: How to decrease the information loss from lag variables?

text: I am using a distributed lag model to analyze a time series data. The duration of study period is 18 years, and the observation is yearly data. When including a 1-year lag effect, the first year of the lag variable becomes missing. Then, a 2-year lag effect makes the first two data of lag variable missing, and so on. I am going to analyze five lag effects in my studies, but five lag variables caused 5 missing data. I assume the multiple imputation may help me overcome the loss of information in these lag variables, but the imputation result is not reasonable. Is there any better idea to impute the missing data in lag variables?


title: Data Acquisition in R

text: R is good at processing data. It is an analytic turbo-ginsu. So LabVIEW is great at getting data. It plugs into anything with electrodes. MatLab has simulink. There are toolboxes for data acquisition. What does "R" have? Does it connect to hardware, or is that too "laboratory" for it? What is the "R" interface to hardware for gathering the data that it so beautifully analyzes? There is ephemeris associated with the method of collection and with the method of analysis that can be mutually informative. If I have video input, and I want to do 2d image tracking, why would I use R? I wouldn't because it can't inform the end-product. It has great libraries that traverse the breadth of the field of mathematics, that are made by the best in the field so they tend to be substantially correct. That doesn't mean anything in a product- centric world. You have to be able to productize your work. MATLAB has the MATLAB compiler. LabVIEW is compiled and has standalone runtimes. Both connect to hardware on one side and the user on the other with good math in the middle. R doesn't have that - or does it? Package list, documentation, and examples are requested for this answer.


title: How to analyze two-way within subjects design with binary data

text: I am interested in testing how users perform on three different user interfaces. On each of the three interfaces, a person will try 12 Tasks. So each participant will do Task 1 three times, Task 2 three times, etc. So, in essence the experiment looks like a 3 x 12 (Interface x Task) within subjects design (where each person does 36 tasks). The sample size will be on the small side (roughly a dozen participants). I am wondering how to analyze completion data, which is binary. If a participant is successful completing a task, then they get a 1. For failure, they score a 0. Ideally, I would like an analysis that is something like an ANOVA so I could see if there are differences in Interfaces, Tasks, or an interaction. If I had to pick though, I am most interested in differences between Interfaces.


title: Using non-random sample to make population estimates

text: The health screening for immigration to the US includes a mandatory chest x-ray screening for all adults, with the aim of detecting signs of tuberculosis. If the x-ray indicates the possibility of TB, the applicant then submits sputum samples for culturing. That's a lot of x-rays - and many of them don't result in positive cultures, especially among older applicants (whose lungs can be funky for many reasons). In a study, 1,500 applicants (out of 24,000 during the study period) were given a new TB diagnostic in addition to the standard screening process.* People with chest x-rays indicative of TB were oversampled to ensure a fair number of applicants with culture-positive sputa. Now, the researchers would like to know (for example): if we had performed this additional diagnostic on all applicants, how many people would have tested positive? _Of those_ , how many would have had a TB-indicative x-ray? _Of those_ , how many would have had positive cultures? Demographic, x-ray, and culture results are known for all applicants; new-diagnostic result is only known for study participants. Two approaches for this study have been proposed.   1. Build a logistic regression on the sample and use to predict for all applicants. There are concerns that the standard errors will be too small, as age, x-ray result, and culture result are quite correlated.   2. Pseudo-bootstrap. Sample 24,000 results from the participants' results (with replacement), assigning a sampling probability to each participant based on the number of all applicants with matching characteristics (demographics, x-ray result, culture result). From the resample, calculate desired quantities. Repeat, calculate confidence intervals for the estimates. Which of these two would you use, if either? Why? I'm happy to clarify any questions. Thanks for your help! *Some study details changed to protect the word count.


title: Poisson regression on binary outcomes with `glmnet`?

text: Because my event is not particularly rare, I am more interested in risk ratios than odds ratios for my current generalized linear modeling problem, which calls for regularized regression. I'm doing regularized regression with `glmnet`. I know that you can obtain risk ratio effects instead of odds ratio effects by fitting a Poisson regression with robust error variance, which is preferred because not using robust error variance inflates your confidence intervals. Well `glmnet` does not provide confidence intervals, so I'm only looking at point estimates. My question is whether I could fit a Poisson regression to binary data using `glmnet` without having coefficients being shrunken to zero that shouldn't because of the way `glmnet` treats error variance. Thank you.


title: Non-significant principal components as predictors in a regression

text: I've a question regarding principal component regression (regression of a DV on principal components). I have 4 components in my PCR and the third component is non-significant as a predictor. What does that mean? It's curious that the "weaker" fourth component however is significant. Can anyone explain this to me?


title: Which parameter on x axis for linear regression?

text: Here's a problem. I'm a geologist, so if you think you can help me... well, go slowly please. I want to predict a continuous variable like porosity from remote sensing data. Let's say I have a measure of the reflectivity of the surface of the earth, densely sampled over an area. And I have rock porosity at some locations. So I'd like to calibrate the former to the latter, and transform reflectivity to porosity. Since porosity _p_ is the unknown, I instinctively plot it on the _y_ -axis, against reflectivity _r_ on the _x_. Then do a linear regression (let's pretend life is this simple for a mo), and get my equation $p = a r + b$, and go off on my merry way. But this feels wrong. Reflectivity depends on porosity (in my hypothesis, and let's say that we know it does in nature). So in fact I should plot reflectivity on the y-axis, vs porosity on the _x_. Then get my linear equation, then rearrange to solve for $p = (r - b)/a$. These approaches give different answers. I found this a bit surprising. Most people I've asked about this say, 'does it matter?'. It does, and the difference can be quite big (several percent). **Which one is the correct approach? Why?** (For what it's worth, my answer is the latter approach, because it 'respects' nature, to anthropomorphize a bit... Not sure I can really explain myself any better). Thanks for any insight anyone can offer,   Matt


title: How to quantify the Relative Variable Importance in Logistic Regression in terms of p?

text: [I already searched the site and I could not find the answer to this question] Given a simple example: A logistic regression model is used to predict whether an online shopper will purchase a product (outcome: purchase), after he clicked a set of online adverts (predictors: Ad1, Ad2, and Ad3). The outcome is a binary variable: 1 (purchased) or 0 (not purcahsed). The predictors are also binary variables: 1 (clicked) or 0 (not clicked). So all variables are on the same scale. If the resulting coefficients of Ad1, Ad2, and Ad3 are 0.1, 0.2, and 03, we can conclude that Ad3 is more important than Ad2, and Ad2 is more important than Ad1. Furthermore, since all variables are on the same scale, the standardized and un-standardized coefficients should be same, and we can further conclude that Ad2 is twice important than Ad1 in terms of its influence on the logit (log-odds) level. But in practice we care more about how to compare and interpret the relative importance of the variables in terms of p(probability of the purchase) level, not the logit(log-odds). Thus the question is: Is there any approach to quantify the relative importance of these variables in terms of p? Thanks a lot!


title: Ill-labeled survey repairable by analysis?

text: Suppose a survey has five items with a six-point Likert scale. It has been rolled out with ill-formed labels: *the (2) was labelled by "somewhat disagree" instead of "disagree" - and vice versa the (3) was labelled by "disagree" *analogous ill-forming for (4) and (5): "somewhat agree" and "somewhat disagree" were mixed up Taken literally, the respondents face a 132546-situation instead of a 123456-situation for _some_ of the items. The begin of the survey was labelled right. Can this be repaired by a clever and statistically sound repair-procedure as part of the data analysis, if the answered questionnaires should form a basis for structural equation modelling? From my point of view: clearly no, the questionnaire must be sent out once again in a sound shape. Even if there is a general conditioning of answering this questionnaire as of the labels were in the right order, "conditioning" is not a 100% concept. So, there are two things we would not know: **how many** people took the questionnaire literally and **who** took it literally. This introduces a "second noise" to the structural equation model, which uses p-values etc. to quantify "noise". Now, "noise" and "second" noise cannot be separated anymore. Even if we knew the relative amount of conditioned answering, we still did not know, who answered in an unconditioned fashion - ending up with the non-separable "second noise" again. Normally I am really not that harsh. Thatswhy: any objections to my point of view - is my position too harsh?


title: Weighted Least Squares with Standardized Coefficients

text: I want to understand how weighted least squares regressions work to implement it in a more complex context. I think I'm a good step into that process, but I'm still wondering what the correct way to calculate standardized coefficients for a weighted regression is. I generated a simple Model with Y as dependent and X1-X3 as independent variables and weights W.               library(MASS)          correlation.matrix <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1),nrow=4)      sample.size <- 500      mean.vector <- c(1,2,3,0)           continuous.data <- as.matrix(mvrnorm(sample.size,mean.vector,correlation.matrix, empirical=TRUE))          X <- continuous.data[,1:3]     E <- continuous.data[,4]     Y <- 2*X[,1]+1*X[,2]+4*X[,3]+sqrt(2)*E+2          W <- runif(500,0,1)     Wd<-diag(W)          Xi <- cbind(rep(1,500),X)      Then I tested different calculation strategies and compared them to results from SPSS. Why SPSS? Because it provides standardized and unstandardized estimates for weighted regressions and the people I'm working for are usually using SPSS. Unweighted regressions (1 & 2) are of course no problem and always equal between my R calculations and SPSS. The unstandardized weighted regression results are also equal between calculation approaches and clear to me (3). However, with (4) I don't quite get why SPSS uses standardized variables as input, weights them and standardize them again (4c). I thought that this is not correct and would have expected that 4b is correct. Any ideas on this issue? 1) Unstandardized regression without and with intercept // Equal to SPSS               solve(t(X)%*%X)%*%(t(X)%*%Y)     lm(Y~X[,1]+X[,2]+X[,3],weights=W)     solve(t(Xi)%*%Xi)%*%(t(Xi)%*%Y)     lm(Y~0+X[,1]+X[,2]+X[,3],weights=W)      2) Standardized regression // Equal to SPSS               Xs <- scale(X)     Ys <- scale(Y)     solve(t(Xs)%*%Xs)%*%(t(Xs)%*%Ys)     solve(cor(X))%*%cor(X,Y)     solve(var(Xs))%*%cov(Xs,Ys)     lm(Ys~Xs[,1]+Xs[,2]+Xs[,3])      3) Unstandardized weighted regression without and with intercept // Equal to SPSS               Xw <- X*sqrt(W)     Yw <- Y*sqrt(W)     solve(t(Xw)%*%Xw)%*%(t(Xw)%*%Yw)     solve(t(X)%*%Wd%*%X)%*%(t(X)%*%Wd%*%Y)     lm(Ys~Xs[,1]+Xs[,2]+Xs[,3],weights=W)      4a) standardized weighted regresion (standardized first, then weighted) // NOT!!! equal to SPSS               Xsw <- Xs*sqrt(W)     Ysw <- Ys*sqrt(W)     solve(t(Xsw)%*%Xsw)%*%(t(Xsw)%*%Ysw)     solve(t(Xs)%*%Wd%*%Xs)%*%(t(Xs)%*%Wd%*%Ys)     solve(var(Xsw))%*%cov(Xsw,Ysw)      4b) standardized weighted regresion (weighted first, then standardized) // NOT!!! equal to SPSS               Xws <- scale(Xw)     Yws <- scale(Yw)     solve(t(Xws)%*%Xws)%*%(t(Xws)%*%Yws)     solve(cor(Xw))%*%cor(Xw,Yw)      4c) standardized weighted regresion (standardized first, then weighted, then again standardized) // Equal to SPSS with Intercept!               Xsws <- scale(Xsw)     Ysws <- scale(Ysw)     solve(t(Xsws)%*%Xsws)%*%(t(Xsws)%*%Ysws)     solve(cor(Xsw))%*%cor(Xsw,Ysw)      Thanks Jan-Michael


title: Name for combining related graphics

text: Often we combine related graphics to help us better understand a problem. We do this for poster talks in research (though we include text). Many info- graphics also do this as do this as do dashboards. The problem is I don't know the proper term for **Combining related graphics** , including static graphs. lattice, faceting, or trellis graphics (A subset of what Tufte called multiples) in which data is split into two or more grouping variables, are not what I'm talking about as each facet is the same type of graph as the others. What I mean is the act of combining related graphs usually without text as we do in a poster. Here's an example in `ggplot2` where someone has combined multiple graphics around a particular event: ![enter image description here](http://i.stack.imgur.com/vebLO.png) I accomplish this with `arrangeGrob`, `grid.arrange` or `par(mfrow=c())` \+ `layout` (sometimes) in R. I'm inclined to call this a dashboard as an info-graphic could be this type of display but an info-graphic may only have a single graph (I tend to think most all graphs are info-graphs). This web site says that Stephen Few says this about dashboards: > A dashboard is a visual display of the most important information needed to > achieve one or more objectives; consolidated and arranged on a single screen > so the information can be monitored at a glance. This seems closest of the terms I'm throwing around but am uncertain. It implies that the display should be on a computer. My description is not limited to computer, for example I could display this on printed paper in a magazine. Perhaps there's a better term(s) than dashboard. **What is the term(s) used for combining different type (i.e., not faceted), yet related, graphs into a single display?**


title: AR process with a constant

text: I am having trouble understanding the estimation of an AR process. In some textbooks, the AR(1) process is defined as follows: $y_{t}=\theta y_{t-1}+ϵ_t$ (which does not contain a constant). So the OLS estimator is biased. I am confused about the cause of the bias. It is explained that $y_{t-1}$ is dependent on $ϵ_{t-1}$ although it is independent of $ϵ_t$. However in linear regression, if the equation does not contain a constant, we cannot make sure the expectation of disturbance is zero. So the OLS estimator is bias without a constant. Does it mean that the OLS estimator is unbiased if I add a constant in AR process?


title: Significance of average correlation coefficients

text: I am using fMRI data (n=25) to predict performance on a behavioral task. At each of 12 time points, I train a linear model on data from n-1 subjects and use this model to predict the left-out subject's performance during a different period of time. For example, I use neural data at time point 1 to predict performance during time window A, and data at time point 2 to predict performance during time window B, etc. For each time point, after leaving each subject out once, I perform Spearman's rank correlation on the observed and predicted values. This gives me a correlation coefficient at each of my time points. I want to report the average prediction across all of my time points, so I Fisher transformed these r values, averaged the resulting z scores, then transformed the average back to an r value. Is this okay to do since my correlations were all computed within the same group of subjects? Based on previous questions I _think_ this is right but please correct me if I'm wrong. My next question is about how to report the significance of this average r. Is this possible? This answer leads me to believe that reporting a confidence interval would be better, but I'm not exactly sure how to do this, or if it's okay because, again, my correlations were computed from data collected at different times from the same group of subjects. Thank you for any advice you may have!


title: studentized residuals in nonlinear regression

text: I am using the nls function in R to perform a nonlinear regression and need to calculate the studentized residuals. Is this something I need to manually do? It seems like I need to manually do a QR decomposition and calculate the hat matrix to do this. I don't have a problem doing this but I wanted to ask if anyone knew an easier way to do this.


title: Feature Selection ARIMAX

text: I am trying to implement the ARIMAX model. I have a lot of exogenous features and I have no clue which ones are important in predicting my endogenous variable. Does ARIMAX have some sort of native feature selection in that it just places less weight to the variables that aren't that important and more weight to the ones that are? If not, should I be using some sort of feature selection algorithm to determine which variables are important?