title: Interpretation of two indexes Interaction Term

text: Respected Fellows. I will thankful if someone help me to explain my model results.my model is as follows. Yit=αPFit+βPSit+δ (PF*PS) it+εit Where Y is GDP per capita PF=Political Freedom Index ranges from (0-10) PS is political stability Index ranges from (0-10) PF*PS=interaction term I=country and t represent time period. I know how to interpret the results of two continuous variables but I am confused how to interpret the results of two indexes interaction term.(for purpose of easiness i exclude other control variables here) My estimated model results as follows. Y= -0.009+0.009PF+0.015PS-0.001(PF*PS) Regard


title: Cox regression when reference group had zero events

text: I would appreciate some advise on an a problem I ran into. I use SPSS for statistical analysis of a study. The study look into how a blood test predicts mortality with patients followup of 1 year. Patients were divided into quartiles, and using the first quartile as a reference group I used Cox regression to determine the unadjusted Hazard ratios (HR) and the adjusted HR for quartile 2-4. The problem is that the reference group (Quartile 1–the first 25%) did not have any events. Therefore I am getting Hazard ratios of 80,000+ for the 2nd, 3rd and 4th quartiles as well as error messages under SPSS. Can anyone advise me on what Im doing wrong?


title: Confused with this confidence interval statistics

text: The formula used to compute a confidence interval for the mean of a normal population when n is small is the following. ![http://i.imgur.com/vluLbC3.gif](http://i.stack.imgur.com/Rm63b.gif) What is the appropriate t critical value for each of the following confidence levels and sample sizes? (Round the answers to two decimal places.) (a) 90% confidence, n = 17 (b) 90% confidence, n = 13 (c) 99% confidence, n = 24 (d) 90% confidence, n = 25 (e) 95% confidence, n = 13 thank you so much in advance!


title: It it possible in R to specify a regression formula for the hazard rate for a survival analysis model?

text: I am currently trying to fit a survival analysis model which has the following survival function: $S(t) = \lambda_i e^{-\lambda_i t}$ but with $\lambda_i = e^{\beta_0 +\beta_1 log(1+X_i)}$ where $X_i$ represents each unique observation. I am trying to use the default survival package but am not sure where I can program in the regression formula for the hazard rate. I believe I might be misunderstanding what the survival analysis package does but would anyone be able to tell me if this can be even done in the R package? Thanks!


title: Using two hazard ratios to make third comparison

text: I have two sets of hazard ratios (with confidence intervals), one comparing risk of developing a cancer between current smokers and never smokers, and another comparing former smokers to never smokers. I do not have the raw data. Is there a way to use this information to create a hazard ratio comparing current smokers and former smokers? For example, the hazard ratio for current smokers vs never smokers is 2.0 [95% CI, 1.2-3.4], and the ratio between former smokers and never smokers is 1.5 [95% CI, 0.9-2.5]. Is it possible to calculate a ratio for current vs former, with a 95% CI? Thanks!


title: What do vertical bars mean in statistical distributions?

text: What do the **vertical bars** mean in the first and third formulae? $$v_i|z_i=k,\mu_k\sim\mathcal{N}(\mu_k, \sigma^2)$$ $$P(z_i=k)=\pi_k$$ $$\pi|\alpha\sim \text{Dir}(\alpha/K1_K)$$ $$\mu_k\sim H(\lambda)$$ This formula is originally from here.


title: Which Test to Use?

text: I was wondering if someone could help me regarding the statistical tests I should use. Basically I have the following: 1) I carried out a survey asking teachers and students to rank components in my conceptual framework using the the scale: (not important) 1 2 3 4 5 (essential) 2) I want to first show what the teachers thought were important and not so important components. My research suggests a benchmark value of 80% is a good benchmark so I want to check if the teachers average response are >= 4 (80% of 5). 3) Next I want to show what the students thought by doing similar tests again checking if average response >= 4 4) Next, I want to compare teachers response to students responses where their responses are different (i.e. where one group thought a component is important whilst the other thought it is not) 5) Finally, I want to be able to decide which components are important and should be kept in the framework and which should be removed based on the responses of the teachers and students. I would be grateful if someone could help me understand what test I should use for each of the above steps. I already did the following: For 2) I did a one-sample t-test with a mu of 4 For 3) I did a one-sample t-test with a mu of 4 For 4) I did a two-sample t-test But I think the above is incorrect since my research indicates:   * that the population must be normally-distributed. Lilliefors (Kolmogorov-Smirnov) normality test suggests it is not.   * Should be interval data. I have an ordinal scale and there are a lot of cases where t-tests are used with oridnal/likert scales.   * The results showed that teachers thought component A was not very important (mu <4), students thought component A was at least very important (mu > 4) and the two-sample t-test showed not significant between the groups. So not sure really what conclusion or further test could be done here. I think I might need to use a Mann-Whitney test but I am not sure so any advice would be greatly appreciated.


title: Weak convergence of an empirical process. Demonstration

text: The empirical process $B_n(x) = \sqrt{n}(F_n(x) − F(x))$ converges weakly to a zero-mean Gaussian process, $B$, with covariance function: $\mbox{cov}(B(x), B(y)) = F(\min\{x, y\}) − F(x)F(y)$. How I can prove this assumption (about covariance function)?


title: Proof/Derivation of Residual Sum of Squares (Based on Introduction to Statistical Learning)

text: On page 19 of the textbook Introduction to Statistical Learning (by James, Witten, Hastie and Tibshirani--it is freely downloadable on the web, and very good), the following is stated: > Consider a given estimate $$\hat{Y} = \hat{f}(x)$$ Assume for a moment that > both $$\hat{f}, X$$ are fixed. Then, it is easy to show that: > > $$\mathrm{E}(Y - \hat{Y})^2 = \mathrm{E}[f(X) + \epsilon - \hat{f}(X)]^2$$ > $$ = [f(X) - \hat{f}(X)]^2 + \mathrm{Var}(\epsilon)$$ It is further explained that the first term represents the reducible error, and the second term represents the irreducible error. I am not fully understanding how the authors arrive at this answer. I worked through the calculations as follows: $$\mathrm{E}(Y - \hat{Y})^2 = \mathrm{E}[f(X) + \epsilon - \hat{f}(X)]^2$$ This simplifies to $[f(X) - \hat{f}(X) + \mathrm{E}[\epsilon]]^2 = [f(X) - \hat{f}(X)]^2$ assuming that $\mathrm{E}[\epsilon] = 0$. Where is the $\mathrm{Var}(x)$ indicated in the text coming from? Any suggestions would be greatly appreciated.


title: What to do when the standard error equals 0

text: I am working on a project where I am evaluating software packages. I have got the results of two experiments where only one variable has been changed. Each experiment has been conducted 10 times. I am doing paired-sample t-tests. For some of my metrics, the results are exactly the same for the ten iterations. Therefore I have a standard error of zero and I cannot do the t-test. If I am correct you can use the p-value generated during a t-test to support your claim that the different values of the variables had an impact on the results (as you can say that your results are statically significant). In my case, what can I say/do to support such claim? Thanks EDIT: From the comment, I am going to add details on my design. I'm am evaluating file carvers, these are digital forensics tools which can be used to extract files from forensics images (such as hard drives). I've planted a known number of evidence (500) and configured each file carver to look for these evidence. I have got tow set of evidence. One is a "well known" evidence (that I created myself for the purpose of the research, so I am planting 500 time the same file) and a realistic set of evidence (composed of 500 unique files). I am doing two tests with the same tool: one for each set of planted evidence. If I understood correctly the different type of t-tests, I can do a paired t-test? From the results that I have got, the tools fail each time to recover all evidence. However, the recovered evidence are the same each time. The tools are working as a search function. Such function might fail if it is run hundreds of times, maybe I need much more. I have tested four tools with different scenarios and the results never varied (same true positives, false positives and false negatives) What I need, is a way to support my findings. When you look at the literature you get told how to talk about p-values results. I just have no idea what I should say if it cannot be computed as their is no variation in my results. I found some articles where the researchers had similar results but they were just ignored (only displayed in the table). It might be obvious but as it is the first time that I'm using stats in my work I'm a bit lost (and I'm self- taught so I might have poorly understood some concepts) Finally, my main question is: The results show clearly that the type of evidence as an impact on the results. How can I support such claim in a scientific fashion?


title: Can I multiply and round proportion data for analysis in a GL(M)M with Poisson errors?

text: Maybe the dumbiest question ever posted on CV: I want to analyze the relationship between proportion data and some covariates in a generalized linear (mixed, but I don't think it matters) model. The natural distribution to use would be binomial but most of the proportions are close to 0 and their distribution looks like a (slightly overdispersed) Poisson. My proportions are by definition not integer values, and so cannot be directly modeled using some implementations of gl(m)m with Poisson errors. For instance I'm using MCMCglmm in R which (rightfully) only allows integer values (and Beta-binomial is not implemented). Could I simply multiply my proportion by let's say 1000 (to keep some accuracy), round them (to make them integer) and do the glmm with Poisson errors on them. Then I could rescale the parameter estimates (and in glmm the _relative_ importance of variance components will stay the same)? I assume (maybe wrongly) that there would be not problem with this, but as I've been unable to find it discussed somewhere, I'm starting to doubt, and would be pleased if others thinking about this later can find an answer (whatever it is) here. Many thanks for the support.


title: Weights in the fit.variogram method of gstat

text: I'd like to know if I understood correctly the following. In the `fit.variogram` method of the gstat library, there is a `fit.method` argument. In the documentation, it says that: > The default method uses weights $N_h/h^2$ with $N_h$ the number of point > pairs and $h$ the distance. This criterion is not supported by theory, but > by practice. For other values of fit.method, see table 4.2 in the gstat > manual. Let's say we have a variogram that includes sets of points separated by some distances $h: 100, 300, 600$ and $900$ and each $N_{h}$ contains 100, 1000, 1200, 900 points, respectively. Does the weighting method described above mean that the weights are calculated as $100/100, 1000/300, 1200/600, 900/900$. That is, we are placing less weight on points that are separated by larger distances. Is this correct? I know it says that this method is not supported by theory, but it seems a bit counterintuitive.


title: Visualization of an ensemble of sequences

text: What are the standard techniques of plotting a statistical ensemble of sequences? I know the chaos game representation:   * H. Joel Jeffrey (1992), Chaos game visualization of sequences. Computers & Graphics 16 (1): 25-33. and two-dimensional plots:   * H.C. Lee Hao and Shu-yu Zhang (2000). Fractals related to long dna sequences and complete genomes. Chaos, Solitons & Fractals, 11:825, free pdf.   * J. Rodriguez-Laguna, PM, M. Ibáñez Berganza, M. Lewenstein, G. Sierra (2011). Qubism: self-similar visualization of many-body wavefunctions, arXiv:1112.3560 [quant-ph]. Are there any others? Edit: I am the most interested in finite alphabets (say, of the order of 2 to 9), with no internal structure and with sequences of any length (also infinite; but in practice 8-18 letters are fine for me). However, other approaches are welcome as well.


title: Is it allowed to transform x and y in an ANCOVA?

text: in my ANCOVA analysis, I have the case that my Y and X show a nonlinear pattern. Transforming either only Y or X does not help to stabilize variance. I did not find any hint that doing an ANCOVA in log-log space is illegal, however, I also did not find any note on that one should try out a transformation of Y AND X... Could someone give me a hint if I am doing wrong with log-log? thanks!


title: Levin Lin Chiu test in stata

text: I need to do a unit root test with my panel data. I have 1500 firms from 98-2000. I want to examine if my LOGwage variable has unit roots or not. When I xtset my panel Firm and year, stata said that it's strongly balanced. But when I do the test xtunitroot llc wage if FirmID, it gives me "Levin-Lin-Chiu test requires strongly balanced data"?


title: Finding the MLE

text: Let $x_1,\ldots,x_n$ be a random sample with pdf $$f(x)= \begin{cases} (\alpha+1)x^\alpha, & 0 \le x \le 1, \\ 0,& \mbox{otherwise}. \end{cases}$$ Find the MLE of $\alpha.$ So, I've found the log-likelihood and taken the derivative, but I'm not quite sure if I've done it right. This is what I have so far. Would I then just set the equation equal to zero and solve for $\alpha$? \begin{align*} \log L &= \sum_{i=1}^n \log(\alpha + 1) + \sum_{i=1}^n \log(x_i^\alpha) \\ &= n\log(\alpha+1) + \sum_{i=1}^n \log(x_i^\alpha)\\ &= \frac{n}{\alpha+1} + \alpha\sum_{i=1}^n \log(x_i) \end{align*}


title: Expected value of multiplication of Identically distibuted random variables

text: i am trying to understand if the following statement is true: $$ E(XY)= E(X^2)=E(Y^2) $$ if $X$ and $Y$ are identically distributed but not necessarily independent r.v. This means that if the variables are identically distributed calculate the expectation of the multiplication of any pair of them is the same that calculate the expectation of one of them squared, due to the fact that they can be treated as _the same_ variable because we are dealing with the expectation of them. Thanks


title: How can I perform a chi-square test to do feature selection in R

text: Document/Term Freq  C1  C2  C3  C4  C5  C6  C7  C8     C9   Tag     D1                  0   1   3   1   1   0   0   1       0   X1     D2                  1   1   3   0   1   0   0   2       0   X2     D3                  2   0   2   0   1   0   0   0       0   X3     D4                  0   0   4   0   1   0   0   0       0   X1     D5                  0   0   1   1   1   0   0   1       1   X2     D6                  0   0   0   0   1   0   0   1       1   X2     D7                  0   0   0   1   1   1   1   3       0   X3     D8                  1   0   0   0   1   2   1   0       0   X1     D9                  1   0   1   1   2   2   1   3       1   X1      I have a DocumentTermMatrix of the above form where D1, D2, D3.....D9 are the documents. C1, C2, C3....C9 are the terms available in all or any of the documents and the numbers are the frequencies that the term occur in the document. The values in the Tag column are the label for each document. My original DTM has around 100000 rows and 2000 features(column) and the number of unique label is around 6000. Now I want to perform chi- square test for feature selection in R. But I don't know how to do that in R and once done how to interpret the results so that I can do a feature selection based on the chi-square result.


title: Moment generating function of a distribution

text: I want to find the moment generating function (mfg) and mean deviation of this distribution: $$f(x,\epsilon,k,\theta) = k\theta^{(1+1/k+\epsilon/k)}x^{(k+\epsilon)}\exp{(-\theta x^k )}/(\Gamma(1+(1+\epsilon))/k)$$ where $\epsilon, k, \theta$ are the three parameters of this distribution. HERE x is from 0 to infinity.


title: Stata equivalent of array

text: I've been trying to solve this problem reading different material but I just can't figure it out. It's pretty simple and takes about 5 seconds to do in languages that work with arrays, but I can't work it out in Stata. Please ignore the substitution of the backticks in some places, I couldn't get them to appear correctly. Let's say I have an ID variable and, among others, variables that show me whether a person has company or is alone and whether he drinks alcohol or not. If someone has company, I have a variable showing me how many these companions are and who they are (by their IDs in the data set). I'm trying to see if these companions drink alcohol and do a repetitive process with that info. Obviously in Stata if the companion has ID 345 I can't use the form alcohol[345] because that gives me the alcohol state of the observation 345, not of the ID 345. Is there a similar form in Stata? I can only list it by doing list alcohol if id==345. What I want to do is take these IDs and pass them on into a repetitive process that allows me to see whether all of these companions don't drink alcohol, part of this process is to have a variable d which I modify for every companion who does not drink alcohol (i.e. when the variable alcohol is not 1). I do the following, if m is the number of companions (let's say two here) and comp'k'`x' are the IDs of the companions of the person with ID k (let's say ID 8 here) :               gen d = .       local m = 2       local k = 8       forvalues x= 1/'m' {         replace d = 'x' if alcohol[comp'k'`x'] !=1       }        so for example the first companion of person with ID 8 would have ID comp81= 345, the second comp82=567. But of course this code does not do what I want, for the reason above. The following will give me the alcohol state of the person whose companions are IDs comp81 and comp82, and not what I want which is the alcohol state of these companions:               . local m=2       . local k = 8       . forvalues x= 1/'m' {       2. list alcohol comp'k'`x'       3. }        And the following doesn't work either:               local m = 2       local k = 8       forvalues x= 1/'m' {       2. list alcohol if id==comp'k'`x'       3. }        I don't just want to list the alcohol state of the companions, but the above are examples to show that I can't find the Stata equivalent of an array and use it to solve the problem I have.


title: What is meant by effective parameters in machine learning

text: My question might be a bit ambiguous, but I started to wonder what does the "effective parameters" mean in machine learning? I have heard few professors of machine learning in my university talk about effective parameters (the context was about k-nearest neighbors or Gaussian mixture models etc.). Any clarification on what this might mean? Thank you for any help! P.S. I have also seen this term in the book: The Elements of Statistical Learning 2nd edition, Trevor Hastie, page 15.


title: Probability function of number of fails to first success

text: I quote the question: > In some cultures, it is important to have at least one son. The plan is to > have children _until one son_ is born. Find the PDF of the number of > daughters in a family. The success rate of having a son is _p_. My first approach: Bernoulli experiment let D = number of daughters. Then d+1 experiments with d successes have been performed. $f_D(d) = {d +1 \choose d} (1-p)^d p^1 = (d+1) (1-p)^d p$ This is apparently **wrong**. The answer should be $f_D(d) = (1-p)^d p$ Why is my approach not correct and what is a correct approach then?


title: complicated model without a huge dataset

text: I've got a set of data that I'm trying to model. Lots of the data is missing, so I'm using multiple imputation. I've got about 360 observations and 13 variables. I'm also using GAMs, but that should be really relevant (I specify the knots in advance). In addition, I've got 4 interactions -- 2 two-way interactions and 1 3-way. All of them are pre-specified by theory, and I see which ones I can cut from the model (using AIC) to achieve more parsimony, after fitting. Here is my problem: I have way too many interactions and variables for the size of my dataset, yet none of them want to come out of the model -- removing them increases AIC. When I compare the model with interactions to the model without, an F-tests suggest that the complicated model is quite significantly different from the simple one. Yet because of high co-linearity, terms jump around quite a bit when I set the seed differently or come across a new data point to add. My prediction results still have the same general gist, but inference on model terms isn't stable. I'm not sure what to do here. Rules of thumb suggest simple specifications for small datasets. But simple specifications are clearly inferior. Beyond "pick a better problem with better data," what would you suggest to someone in this situation? My gut instinct is just to de-emphasaze inferential results and stick with prediction when I write this up. But that feels like a coverup. Suggestions welcome. Apologies for the vague question.


title: Birthday problem, but matching couples instead of individuals

text: Suppose _n_ couples are invited to a party. What is the probability that there are at least two husband–wife pairs such that the husbands have the same birthdays and so do their wives?


title: Using a normal distirbution for count data

text: I was just wondering if it is possible to use a normal mixed model for count data. Would it be better to use a Poisson regression if your response is count data?


title: How to account for aggregation over subjects in meta-analysis?

text: I conduct a small meta-analysis of the preferences ratings for objects in two different conditions. I compute Fisher's z and SE for each study. 1) There are studies with totally different number of trials. For example, there is a study with 24 subjects and 2 trials per subject and a study with the same 24 subjects and 44 trials per subject. How to account for such variability? The formula for SE of z (1/sqrt(n-3)) is inversely dependent on a sample size, but what to consider as a "sample" in such situation? The number of subjects, the number of trials, or a multiplication of both, i.e. the total number of trials per experiment? 2) There are studies that use non-aggregated data. For example, in one study a binomial "like vs. don't like" scale is used and a total amount of "like" and "don't like" answers is reported. To compute z I first compute Cox's d, then transform it and its variance to z and SE of z. The resulting SE of z is quite small in comparison to a value computed using 1/sqrt(n-3). How to correctly include this study in analysis?


title: Articles that describe formulas of how to calculate confidence intervals for the omega-squared effect size?

text: I'm looking for articles that describe methods of how to calculate CIs for the omega-squared effect size from an one-way ANOVA. The closest I've come is Finch & French (2011) _A Comparison of Methods for Estimating Confidence Intervals for Omega-Squared Effect Size_ link, but they don't present any formulas. I'd also welcome any links to an R implementation of either analytical estimates or bootstrap CIs. Thanks for any help on this matter!


title: Advice on the statistical analysis of clinical data using summary statistics, t-tests, ANOVA, and linear regression

text: I need some help with the statistical analysis of a study of a particular surgery to remove a particular cancer. I am using the statistical program R to conduct my analysis. My data are saved in the object `study_data`. ### Data               # Create reproducible example data     set.seed(50)          study_data <- data.frame(                   Patient_ID = 1:500,                   Institution = sample(c("New York","San Francisco","Houston","Chicago"),500,T),                   Gender = sample(c("Male","Female"),500,T),                   Race = sample(c("White","Black","Hispanic","Asian"),500,T),                   Tumor_grade = sample(c("One","Two","Three","Four"),500,T),                   Pathologic_stage = sample(c("P0","Pa","Pis","P1","P2a","P2b","P3a","P3b","P4a","P4b"),500,T),                   Treatment_arm = sample(c("One","Two","Three","Four"),500,T),                   Surgery_age = round(runif(500,20,100)),                   Nodes_removed = round(runif(500,1,130)))      Here is what the data look like:               # Peak at the first six lines of the data     head(study_data)            Patient_ID   Institution Gender     Race Tumor_grade Pathologic_stage Treatment_arm Surgery_age Nodes_removed     1          1       Houston   Male Hispanic         One              P2b           Two          77           130     2          2 San Francisco Female Hispanic       Three               Pa           Two          38           112     3          3      New York Female    Black        Four               P0          Four          90            90     4          4       Chicago   Male Hispanic         Two              Pis          Four          46             4     5          5       Houston Female    Black        Four              P2a          Four          96           114     6          6      New York   Male    Black       Three              P3b          Four          92             7      ### My interest I am interested in learning more about what variables are associated with the number of lymph nodes removed during the surgery. My first thought was to simply stratify the data by a particular variable and then calculate the median number of nodes removed. For example, to see if the institution at which the surgery was performed mattered, I could write:               cbind(do.call(rbind, by(study_data$Nodes_removed, study_data$Institution, summary)))                        Min. 1st Qu. Median  Mean 3rd Qu. Max.     Chicago          1   25.50   65.5 64.48   98.75  129     Houston          1   40.00   71.0 69.26  100.00  130     New York         4   36.00   67.0 67.96  100.00  129     San Francisco    3   36.75   61.0 65.76   99.00  127      This lets me compare the median nodes removed in each institutional city. ### My question I would like to fully examine the association between all of my variables and the outcome `Nodes_removed`.   1. Should I just do these simple summary statistics for all of my variables?   2. Do I need to perform some sort of hypothesis test for all of the associations to say whether or not the summary statistics differ? For example, should I calculate a median and a confidence interval for each comparison?   3. Or should I be using t-tests to compare one group to another?   4. In the case of a multi-level variable, should I use ANOVA?   5. Is there any role for linear regression analysis here?    6. If I wanted to build a single model that includes every possible predictor variable, what method should I use? For example, say that I am most interested in the association between the age at which the surgery was performed, `Surgery_age`, and `Nodes_removed`. However, I would like to adjust this association for potential confounders like gender, race, tumor grade, treatment arm, etc. What is the best way for me to do this? Thanks for any advice you can give!


title: How to interpret autocorrelation of residuals and what to do with it?

text: I was wondering what does it mean when time series residuals have autocorrelation? How should I deal with it?


title: Language modeling: why is adding up to 1 so important?

text: (if this venue is inappropriate, feel free to migrate it) In many natural language processing applications such as spelling correction, machine translation and speech recognition, we use language models. Language models are created usually by counting how often sequences of words (n-grams) occur in a large corpus and normalizing the counts to create a probability. To account for unseen n-grams, we use smoothing methods (see several listed here) which take some of the probability mass from the n-grams which are attested in the model and distribute this mass among lower order n-gram (shorter word sequences) backoff probabilities. Many of the smoothing techniques become mathematically complex because of the constraint that the calculations must keep the distribution as a probability (must add up to 1). What is the reason for this constraint? What is the advantage of using strict probabilities for prediction instead of scores of any other kind?


title: Possible range of $R^2$

text: Suppose are three time series, $X_1$, $X_2$ and $Y$ Running ordinary linear regression on $Y$ ~ $X_1$ ($Y = b X_1 + b_0 + \epsilon$ ), we get $R^2 = U$. The ordinary linear regression $Y$ ~ $X_2$ get $R^2 = V$. Assume $U < V$ What's the minimum and maximum possible values of $R^2$ on regression $Y$ ~ $X_1 + X_2$ ($Y = b_1 X_1 + b_2 X_2 + b_0 + \epsilon$ )? I believe the minimum $R^2$ should be $V$ + a small value, since adding new variables always increases $R^2$, but I don't know how to quantify this small value, and I don't know how to obtain the maximum range. Thanks a lot!


title: Can just about any model be estimated using maximum likelihood?

text: If someone wants to estimate a non-standard model -- for instance, one where the functional form is not explicitly spelled out in statistics textbooks -- can most such models be estimated using maximum likelihood optimization? The assumption would be that the residuals are normally distributed, and then the model would be estimated recursively with the conditional log likelihood given, for example, by James D. Hamilton in "Time Series Analysis" on page 132. The standard errors would then be taken from the estimate of the Hessian matrix. I could see a few situations, such as irregular data sets, where this might be tricky. But this is just a more general question.


title: A framework for comparing the performance of Expectation Maximization

text: I have my own implementation of the Expectation Maximization (EM) algorithm based on the following paper http://pdf.aminer.org/000/221/588/fuzzy_k_means_clustering_with_crisp_regions.pdf I would like to compare the performance with another implementation. For the testing, I am using K number of centroids with 1 Gb of txt data and I am just measuring the time it takes to compute the new centroids in 1 iteration. I tried with the EM implementation in R, but I couldn't, since the result is plotted in a graph and get stuck with some large number of txt data. Does anybody know of an implementation of EM to measure the performance?


title: Are SVD (Singular Value Decomposition) values always positive? Is there a relation between the maximum SVD value and the original data?

text: Assuming it's the standard SVD (no variation of it) with $A = USV^T$, would the $A$ matrix always have positive values (0 to $\infty$)? I noticed that the $U$ and $V^T$ matrices had some negative values with the sample data I used, but I want to be sure that the $A$ matrix has only positive values so that I can choose the proper normalization technique. Also, is there is a mathematical relation between values in the original matrix and values in the $A$ matrix of the SVD? For example, if your original data's mean and range are so and so, then the maximum value in the SVD will be some function of that.


title: Finding post-hoc power of a repeated measures ANOVA using G*Power (Am I entering the right parameters?)

text: I would like to report the power observed with my mixed model analysis (or repeated measures ANOVA). My experiment followed a completely within subject design with 2 independent factors (two levels & 5 levels, respectively). 12 participants were recruited to run the 10 combos (5x2) twice (repeated observations), yielding 240 observations. There were some missing data, so the actual total observations was 230. I have never used G*Power before but it seems like a pretty powerful tool. I attempted to run a post-hoc analysis as shown below. The resulting power is so high (0.933) that I'm a little skeptical. I've taken the sample size to mean # of subjects (12) and that the number of groups refers to levels for between- subject factors, which I have none. The corr among rep measures was just something I estimated [edit: I thought this parameter isn't making a big difference, but it is. Changing it to 0.5 actually lowers the power to 0.5] Did I enter anything wrong? I would appreciate it if anyone could help me figure out this power analysis.               Analysis:   Post hoc: Compute achieved power     Input:  Effect size f   =   0.2         α err prob  =   0.05         Total sample size   =   12         Number of groups    =   1         Number of measurements  =   10         Corr among rep measures =   0.8         Nonsphericity correction ε  =   1          Output: Noncentrality parameter λ   =   24.0000000         Critical F  =   1.9758061         Numerator df    =   9.0000000         Denominator df  =   99.0000000         Power (1-β err prob)    =   0.9333319


title: Time spent in an activity as an independent variable

text: I want to include time spent doing something (weeks breastfeeding, for example) as an independent variable in a linear model. However, some observations do not engage in the behavior at all. Coding them as 0 isn't really right, because 0 is qualitatively different from any value >0 (i.e. women who don't breastfeed may be very different from women who do, even those who don't do it for very long). The best I can come up with is a set of dummies that categorizes the time spent, but this is a waste of precious information. Something like zero-inflated Poisson also seems like a possibility, but I can't exactly figure out what that would look like in this context. Does anyone have any suggestions?


title: Which test to use to compare proportions between 3 groups?

text: We are testing an e-mailing marketing campaign. On our initial test, we sent out two different e-mail types and had a third control group that did not receive an e-mail. Now we are getting back "results" as proportion of users who returned to our app. Here are the results:               Group | received e-mail | returned | %-returned     A | 16,895 | 934 | 5.53%     B | 17,530 | 717 | 4.09%     C | 42408 | 1618 | 3.82%      It looks like Group A may actually be better than B and C, but what is the proper test to show this?


title: Variance associated with factors in GLS (nlme)

text: First time posting here, so thank you ahead of time for your help. I'd like to estimate the variances associated with two factors in a relatively simple, but unbalanced GLS model, and I am unsure how to extract that information (or how to manually calculate it). For the two factors I'm using, DS3 and DS4, DS3 has four levels (WW - ZZ) and DS4 has seven levels (AA - GG). I've been using the nlme package; the model itself and model output are below.               > library(nlme)     > temp_gls4 = gls(output~DS3+DS4,data=temp_month,weights=varIdent(form=~1|DS4))     > summary(temp_gls4)     Generalized least squares fit by REML       Model: output ~ DS3 + DS4        Data: temp_month           AIC      BIC    logLik     1199.769 1268.615 -582.8845          Variance function:     Structure: Different standard deviations per stratum     Formula: ~1 | DS4      Parameter estimates:            AA        BB        CC        DD        EE        FF        GG      1.0000000 1.2900663 1.1351448 0.9501666 0.7463385 0.8933227 0.9123444           Coefficients:                     Value Std.Error   t-value p-value     (Intercept)  3.853300 0.1232730 31.258279  0.0000     DS3XX       -0.053908 0.1505280 -0.358129  0.7204     DS3YY       -0.119480 0.1505280 -0.793737  0.4278     DS3ZZ       -0.146255 0.1552625 -0.941988  0.3467     DS4BB        0.934175 0.2081654  4.487655  0.0000     DS4CC       -0.093199 0.1098115 -0.848717  0.3965     DS4DD        0.831305 0.2011411  4.132945  0.0000     DS4EE        1.867642 0.2538358  7.357675  0.0000     DS4FF        1.429805 0.1680183  8.509820  0.0000     DS4GG        0.050992 0.1779757  0.286510  0.7746      I know that in a traditional balanced least squared situation you can pretty easily calculate the variance of DS3 and DS4 by manipulating the mean squares (e.g. var(DS3)=(MS(DS3)-MSE)/7), but the data here aren't balanced and I don't know the math behind the GLS model and REML as compared to an traditional least squared model. I've tried using the anova command, but it just produces the F statistics associated with each factor. Any help you can provide would be appreciated! Jon


title: What are good RMSE values?

text: Suppose I have some dataset. I perform some regression on it. I have a separate test dataset. I test the regression on this set. Find the RMSE on the test data. How should I conclude that my learning algorithm has done well, I mean what properties of the data I should look at to conclude that the RMSE I have got is good for the data?


title: Can I include future index into current forecasting in R models

text: I am new to R. Now my team is building forecasting models for monthly sales. Our sales correlates quite well with a industry index. As a forward index, we can get both past 5 years' index and next 3 years' data. Now we are using auto.arima to generate models. We'd like to use xreg in R to correlate with this index. However we may need to handle this problem: if this index predicts next month the sales will decrease, how to include this future index into arima model and generate different forecasting result? Are there any other models in R can help me to include this factor? Thanks!


title: P-value of a survival ROC c-index

text: It is possible to calculate the c-index for time dependent outcomes (such as disease) using the survivalROC package in R. My question is : is it possible to produce a p-value for the c-index that is calculated (at a specific point in time)? If I were to use a standard ROC I understand that I can use Wilcoxon test to calculate the p-value of the AUC (here's a simple example). However, because this is survival data I don't think I can use quite the same approach. I have thought of just calculating the concordance index using the survcomp package and getting the p-value from that, however `survivalROC` predicts the c-index at a specific time, whereas the concordance index doesn't seem to. Any suggestions?


title: How to, or what is the best way, to apply propensity scores after matching?

text: I have developed a recent interest in propensity scores. I have been using the SPSS tool created by Dr. F. Thoemmes to calculate propensity scores using bivariate "treatment" variables (e.g., depression) and several covariates (e.g., age, sex, persons in household). I am then given a resulting propensity score, but am left wondering what to do with it. I have read that what is typical would be to match two individuals who have nearly identical propensity scores (e.g., two people get group number "17") but who actually differ on your treatment variable (e.g., depression), and then do a paired t-test based on group number and dependent variables (e.g., household income). In this example, we would see how two individuals with all of the same propensities (e.g., age, sex, persons in household) but differing on your treatment variable actually differ on your dependent. This idea makes sense to me, but the software actually does not do matching based on propensity scores, and I don't know how to match them using SPSS or Excel, and I don't want to currently bother to learn how to do so in another program/language (e.g, R). This laziness, lets call it, has forced me to do more research. Two authors state: "After the matching is completed, the matched samples may be compared by an unpaired t-test. (“Matching” erroneously suggests that the resulting data should be analyzed as if they were matched pairs. The treated and untreated samples should be regarded as independent, however, because there is no reason to believe that the outcomes of matched individuals are correlated in any way)." (Schafer & hang, 2008). Other research seems to suggest people often input propensity scores in logistic regressions next to their independent variable of interest, and see how the independent variable predicts while "propensity" is controlled for. Although this line of research is interesting, I must admit I am slightly lost regarding what methods are possible/best in terms of conducting quantitative analyses AFTER propensity score calculations. Any guidance on this matter would be appreciated. I will likely have follow-up questions, too! EDIT: I want to highlight that I am concerned about what types of inferential analyses to do AFTER I get the propensity scores calculated for each individual. For instance, perhaps I could calculate propensity score of being depressed (yes,no) based on covariates (age, number of people in household, smoking, sex, state). The program calculates a propensity score as a new variable for each individual. AFTER this, I am interested in seeing if depression is associated with household income while "controlling for"/"considering"/"matching" (Chose word based on what method you suggest, perhaps) the effect associated with propensity.


title: What is "empirical regression"?

text: I was asked how much I know of "empirical regression." I have never heard the expression. A web search yielded nothing useful. I suspect it is a term coined by someone to refer to some ad hoc procedure which might be known by other name. Has anyone have a reference?


title: Two-Sample t-Test for Equal Means with unequal variances for large samples

text: How can I perform a two sample test of means with unequal variances for a very large sample in R? In case of large samples the statistic will asymptotically follow a normal distribution. Which R function will help me to do this?


title: Logistic regression and complementary log log model

text: is there like a diagnostic checking(like those in the time series) done for logistic regression and complementary loglog model?


title: How to find $E[x|y]$ when distributions of y and x are separately known,(p.s. they are both Gaussian)?

text: In detail, I have these relations (in order of causality): $u_1 = ax_0$ $x_1 = u_1 + x_0$ $y = x_1 + w$ where $w = N(0,1), x_0 = N(0,\sigma^2)$. This was my approach: I know the distribution of $x_1 = N(0,(1+a)^2\sigma^2)$ and I know $y = N(0,(1+a)^2\sigma^2 + 1)$ The expectation to find is: $E[x_1|y]$. What I'm unable to figure out is how to go about computing this expectation, and the inter-dependence of the two variables: $x_1$ and $y$. I know that $y$ depends on $x_1$ but is the converse also true?


title: Combining two separate Cox PH models into one model?

text: I have two Cox proportional hazards models (in R), using same outcomes and predictors, one for $n_m$ males and one for $n_f$ females. Is it possible to combine them into one equivalent model over all $n_m+n_f$ individuals, using for example an interaction between the predictors and the sex? For example, this doesn't work in the sense of having same coefficients as each of the sex-specific models:               library(survival)     n <- 200     d <- data.frame(time=rexp(n), event=sample(0:1, n, replace=TRUE),        x=rnorm(n), sex=sample(c("m", "f"), n, replace=TRUE))     l1 <- coxph(Surv(time, event) ~ x, data=subset(d, sex == "m"))     l2 <- coxph(Surv(time, event) ~ x, data=subset(d, sex == "f"))     l3 <- coxph(Surv(time, event) ~ x * sex, data=d)          > coef(l1)          x     -0.1841861     > coef(l2)          x     0.01391554     > coef(l3)            x         sexm       x:sexm     -0.009145576  0.060432017 -0.179186702


title: Simple cloud computing to run R + JAGS simulations

text: I want to simulate the frequentist properties of a Bayesian model. So, for example, I might want to fit a Bayesian model 1,000 times to 50 different configurations each of which takes about 10 seconds to fit on my machine. i.e., total computing time = 1000 * 50 * 10 / 60 second / 60 minutes / 24 hours = 5.7 days on my machine. Running this on my machine is not particularly practical. I could reduce some of the features of my simulation (e.g., fewer simulations per cell of the design; fewer cells in the design; shorter chains for each simulation). However, it would be nice if I could ship this simulation off to the cloud. I could then run it and collect the results when it was finished. The simulation is also highly parallel, such that it could potentially take advantage of using multiple computers to speed up the time it takes to run. **Questions:**   * How can I get a ballpark on the costs of such analyses? For example, if I see a figure that computing time costs 10 cents per hour on EC2, can I assume that the performance of the EC2 machine would be similar to a modern desktop machine.   * Are there any cloud computing systems particularly suited to running occasional R and JAGS simulation jobs?   * Are there any trade-offs between running the jobs in parallel on perhaps multiple instances versus running the program on just one instance?   * Are there any tutorials or examples of running R and JAGS in the cloud? ### Initial discoveries I found this post that provides some useful information about using JAGS on EC2.


title: using stl (seasonal decomposition with loess) for weekly data

text: I am trying to decompose a weekly time series using the R function 'stl'. One of the important argument of this function is the number of data per cycle. Naturally in this case one would choose 52. However, my data is published every week, on the same weekday. Hence I have sometimes 53 data per year. Not a big deal, I agree. But I would like to be precise here. Is there a commonly accepted methodology to treat these situations? I can think of several solutions: 1/ I average the last 2 data in years with 53 data?   2/ I average the last 2 or the first 2 data in years with 53 data? (depending on when we are in the year)   3/ I add an extra datapoint in years with 52 years?


title: Manually evaluating the DIC: very big number of effective parameters?

text: My problem is the following: I'm evaluating the fit to a function $f(\mathbf{x},\theta)$ via MCMC (because I have some priors on the parameters), and I'm trying to evaluate the DIC, given by: $$\rm{DIC}=\bar{D}+p_D,\ \ \ \ (1)$$ where, if we define the deviance $D(\theta)=-2\log(\mathcal{L(\theta|\mathbf{x})})$, and $L(\theta|\mathbf{x})$ is my likelihood, then $$\bar{D} = E^{\theta}[D] = -2E^{\theta}[\log(\mathcal{L(\theta|\mathbf{x})})],\ \ \ \ (2)$$ where $E^\theta[\cdot]$ represents the expected value over the posterior distribution, and the **efective number of parameters** , $p_D$, is given by $$p_D = \bar{D}-D(\hat{\theta})\ \ \ \ (3)$$ where $\hat{\theta}$ is the posterior expectation of the parameters. This is so far what I understood from the paper by Spiegelhalter et al. (2002). The thing is that, just for testing purposes, I'm using a simple gaussian likelihood of the form: $$\mathcal{L}(\theta|\mathbf{x}) = \frac{1}{(2\pi \sigma^2)^{n/2}}\exp\left(-\sum_{i=1}^n\frac{\mathbf{x}_i-\mathbf{f}_i(\theta)}{2\sigma^2}\right),$$ and I'm getting **HUGE** values for $p_D$. Considering $\sigma$ as a parameter and the fact that $\theta$ is 12-dimensional in my particular application, I would expect values of $p_D$ close to 13, but I get values in the order of 900! What I'm doing to evaluate everything is the following:   1. After thinning my MCMC chain, for each link, I get a value of the likelihood. This gives me a sample of the value of the likelihood at each link, and then I estimate $E^\theta[D]$ as: $$\hat{D} = \frac{1}{L}\sum_{i=1}^LD(\theta_i),$$ where $L$ is the number of links and $\theta_i$ is the value of the parameters (including $\sigma$) at that link.   2. I get the posterior mean of my parameters from all my MCMC links (supposing I'm sampling from the posterior) in order to get $\hat{\theta}$. With estimates for this and for $E^\theta[D]$, I just replace the values in eqs. (2) and (3) to obtain (1). With this procedure I get reasonable values for my estimate of $E^\theta[D]$, but $p_D$ just doesn't seem right. Is there anything wrong in my procedure? Thanks in advance for the help!


title: How can I compare two ecdfs for p-values from Kolmogorov Smirnov tests if the number of samples in the two tests are different?

text: The plot below shows the ecdf of the p-values from to different statistics, let us identify the blue graph as vector x and the red graph as vector y. Now I execute a ks.test for x and y compare these vectors with a uni distribution in a small range.               ks.test(x[x<=0.2],function(x){punif(x,0.01,0.2)},alternative = "two.sided")     ks.test(y[y<=0.2],function(x){punif(x,0.01,0.2)},alternative = "two.sided")      Now I get better values for the blue-line but it's obviously that the red-line should be a better fit. If I watch closer to the values in x and y I get only 26 values in x and more than 1500 in y. So I think, the reason why the ks.test gives a better p-value for the blue-line depends on the low number of relevant points in x. The values in x and y comes from a simulation. Schould I increase the number of simulated experiments? ![ecdf](http://i.stack.imgur.com/QYHp9.png)


title: Averaging correlation values from data with multiple levels?

text: At least two questions (1, 2) have addressed the question of how to handle the analysis of multiple correlations. I would like to extend this question to the case of data that has multiple levels (say as a relatively straightforward case correlations between personality traits X1, X2, and X3 on trait Y as nested within 40 people). How then would average correlations be best calculated (and tested for significance) on the understanding that a) the scales have been appropriately coded such that sign does matter and b) that the level of inference desired is to the typical person? Would the optimal analysis be in the family of intraclass correlations, or is this something else entirely?


title: Repeated measures across time, two treatment groups

text: I'm having difficulty modeling something… I think it's fairly simple, but I'm just not confident enough in my ability. I setup an experiment with two treatment groups, six replicates in each treatment group, and measured "diversity" at three time points. I want to see whether "treatment" has an effect and if so, whether it has an effect across all timepoints. I have labelled an id for each of the replicates so that I can tell the model that I've repeatedly measured that replicate. I want to do a repeated measures model and this is what I think I have to do:               lmer(diversity~treatment*time+(1|id)+(1|time)      Could someone tell me whether this is correct, or tell me how to find the correct model?


title: Notation of orthogonality

text: Consider three random variables $X,Y,Z$, which are jointly normally distributed. I know that $Y$ is orthogonal to $X$ conditional on $Z$, in the sense that $\beta_{YX;Z}=0$ (i.e. the regression coefficient of Y on X conditional on Z is zero). I'd like to know if it is correct to state this orthogonality condition as: $$ X \bot Y | Z $$ In other words, I'm interested in understanding the precedence of the operator for orthogonality.


title: Survival Analysis in R with Grouped Data

text: I'm just getting started with survival analysis and I'm having trouble finding something in `R` that will do what I'm looking for. Most of the packages use survival objects that have an individual record for each patient being looked at (so, patient $1$ at time $0$ was _alive_ , patient $2$ at time $0$ was _alive_ , etc.). My data however is _grouped_ because of the larger numbers. For example, at time $0$, $1000$ patients were _alive_ and $0$ were _dead_ , at time $1$, $975$ patients were _alive_ and $25$ were _dead_.               (d1 <- data.frame(TIME=10:5, DEAD=c(0, 0, 1195, 1237, 1251, 1257),                       ALIVE=c(1398, 1398, 203, 161, 147, 141)))       TIME DEAD ALIVE     1   10    0  1398     2    9    0  1398     3    8 1195   203     4    7 1237   161     5    6 1251   147     6    5 1257   141      Is there a `package` or something in `R` or a specific model or something that allows me to perform a survival analysis on data formatted this way? I know this question may not be specifically statistical in nature, but I figured people here would know the most about `R`.


title: Null GLMM Poisson underestimate the mean of the response variable. Is it indicative of poor fitting?

text: I want to test the fixed and random effects of some covariates on a discrete variable with non negative values. In exploratory analysis I fitted a null Poisson GLM and an null Poisson GLMM. However, the GLMM underestimated the mean value of the response variable even after inclusion of fixed and/or random covariates. I also tried Bayesian approaches, zero-inflated models and negative binomial distributions but the "problem" remains. Response variable mean: 0.7804   GLM intercept: 0.7803772   GLMM intercept: 0.6595108 Is the estimated intercept of the GLMM an indicative of poor fitting of the model?               summary(banco2$caes)          Min. 1st Qu.  Median    Mean 3rd Qu.    Max.         0.0000  0.0000  0.0000  0.7804  1.0000 12.0000                 mod1 <- glm(caes ~ 1, poisson, banco2)       Deviance Residuals:           Min       1Q   Median       3Q      Max         -1.2493  -1.2493  -1.2493   0.2381   6.5689       Coefficients:                   Estimate Std. Error z value Pr(>|z|)           (Intercept) -0.24798    0.01078  -23.01   <2e-16 ***       ---       Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       (Dispersion parameter for poisson family taken to be 1)           Null deviance: 15304  on 11027  degrees of freedom       Residual deviance: 15304  on 11027  degrees of freedom       AIC: 27654       Number of Fisher Scoring iterations: 5            exp(mod1$coefficients[1])       (Intercept)         0.7803772                 (mod2 <- lmer(caes ~ 1 + (1 | setor), poisson, data = banco2))       Generalized linear mixed model fit by the Laplace approximation       Formula: caes ~ 1 + (1 | setor)          Data: banco2          AIC   BIC logLik deviance        13575 13590  -6785    13571       Random effects:        Groups Name        Variance Std.Dev.        setor  (Intercept) 0.39817  0.63101       Number of obs: 11028, groups: setor, 559       Fixed effects:                   Estimate Std. Error z value Pr(>|z|)       (Intercept) -0.41626    0.02937  -14.18   <2e-16 ***            exp(fixef(mod2))       (Intercept)         0.6595108        Best regards!


title: Constrained mean variance optimization

text: I have some numbers $X={x_1,x_2,...,x_n}$. I add some small numbers to them and get $Y={y_1,y_2,...,y_n}$ where $y_i=x_i+\epsilon_i$. How can I find $Z={z_1,z_2,...,z_n}$ to minimize $\sum_{i=1}^n(z_i-y_i)^2$ such that $\hat{\text{mean}}(Z)=\hat{\text{mean}}(X)$ and $\hat{\text{var}}(Z)=\hat{\text{var}}(X)$. I have solved it for the simple case of $\hat{\text{mean}}(X)=0$ and $\hat{\text{var}}(X)=1$. However, the computation is too long for the general case. I wonder if there exist more statistical (hopefully shorter) version of the solution. Regards PS: $x_i\in R$


title: A priori comparisons in a repeated measures design

text: I'm interested in finding out more about the use of _a priori_ planned contrasts within a one-way repeated measures ANOVA design. Predictions:   * Grps 2+3 < Grp 1   * Grp 2 < Grp 3 My questions are:   1. **If the data violates RM ANOVA assumptions (e.g., normality, sphericity), how does this affect the accuracy of the planned contrasts?**   2. Following on from this, if transformations were considered but not successful, **is there any way of performing the planned contrasts within a non-parametric alternative?** My sense is 'no', but where do others suggest to go from here? (e.g., is there a way of somehow combining variables in SPSS (e.g., for the first prediction where there are two groups in one cluster)?


title: Trying to understand an example where unbiased estimators don't exist

text: I am new to statistics especially in the topic of estimators and sufficient statistic. I am reading a note which says " **unbiasedness is a desirable (but not necessary) property of a good estimator** ". Then it gives an example where unbiased estimators don't exist. Since there is no further explanation in the example, I am struggling to understand it. The example says: > If $X\sim B(\varphi(\theta))$ for some function $\varphi$, then $\theta^*\in > K_0\iff\mathbb{E}_\theta\theta^*\equiv\theta^*(0)(1-\varphi(\theta))+\theta^*(1)\varphi(\theta)=\theta,\forall\theta\in\Theta)$, > (where $K_0$ is the class of all unbiased estimators). This is just for reference:   An estimator $\theta_0^*=\theta_0^*(X)$ from a class $K$ of estimators of $\theta$ is called **efficient** in $K$ if, for any $\theta^*\in K$, $$\mathbb{E_\theta}(\theta_0^*-\theta)^2\le\mathbb{E_\theta}(\theta*-\theta)^2,\forall\theta\in\Theta.$$ For a function $b=b(\theta),\theta\in\Theta$, let $$K_b=\{\theta^*:\mathbb{E}_\theta\theta^*=\theta+b(\theta),\forall\theta\in\Theta\}$$ be the class of all estimators with the bias $b(\theta)$. I understand the basic idea of estimators, the property for being biased and unbiased. But I have no clue why unbiased estimators don't exist in the given example above. May anyone kindly explain it to me please? Many thanks!


title: Relation and difference between time series and regression?

text: What are relation and difference between time series and regression? For **models and assumptions** , is it correct that the regression models assume independence between the output variables for different values of the input variable, while the time series model doesn't? What are some other differences? For **methods** , from a website by Darlington > There are a number of approaches to time series analysis, but the two best > known are the regression method and the Box-Jenkins (1976) or ARIMA > (AutoRegressive Integrated Moving Average) method. This document introduces > the regression method. I consider the regression method far superior to > ARIMA for three major reasons I don't quite understand what the "regression method" for time series is on the website, and how it is different from the Box-Jenkins or ARIMA method. I appreciate if someone can give some insights on those questions. Thanks and regards!


title: Likelihood ratio test in R for categorical variables

text: I am working with behavioral data of male sea lions, with a binomial model to understand the effect of different variables in determining the location where the encounters between males occur (Land -vs- water). For this I have several quantitative and categorical variables. I am testing the significance of the odds ratio for each variable using likelihood ratio tests. For my dataset the variable "Temperature" was divided into 3 different categories: "1" values<20ºC, "2" values between 20ºC and 30ºC, and "3" values>30ºC. The results of the odds ratio shows that the odds of an aggressive interaction occurring on the land increase in 9% when going from the first temperature interval (<20ºC) to the second interval (20ºC-30º). However, the odds of an aggressive interaction occurring on land decrease in 8% when going from the second temperature interval (20ºC-30º) to the third interval (>30ºC). When applying the likelihood ratio test (lrtest function) to these data it gives me a single value of P, but what I need is one value for each change of the variable: one for the change between the first and the second interval, and one for the change between the second and the third interval. It means I need the likelihood ratio test to include the categorical variable, to understand the significance of each odds ratio. The script I used is this one:               # MODEL WITH ALL VARIABLES          > J_mod.a <- glm(as.factor(Place) ~ as.factor(TempF) + as.factor (AgType) + FemDen,  data=JUNE, family=binomial("logit"))     > J_mod.a          Call:  glm(formula = as.factor(Place) ~ as.factor(TempF) + as.factor(AgType) +      FemDen, family = binomial("logit"), data = JUNE)          Coefficients:     (Intercept)   as.factor(TempF)2   as.factor(TempF)3  as.factor(AgType)2        1.34820             0.08229            -2.53280            -1.73712        FemDen        0.11443            Degrees of Freedom: 124 Total (i.e. Null);  120 Residual     (7 observations deleted due to missingness)     Null Deviance:      159.6           # MODEL WITHOUT TEMPERATURE          > J_mod.TEMP <- glm(as.factor(Place) ~ as.factor (AgType) + FemDen, data=JUNE,   family=binomial("logit"))     > J_mod.TEMP          Call:  glm(formula = as.factor(Place) ~ as.factor(AgType) + FemDen,      family = binomial("logit"), data = JUNE)          Coefficients:     (Intercept)  as.factor(AgType)2              FemDen        0.01463            -1.77035             0.08361            Degrees of Freedom: 124 Total (i.e. Null);  122 Residual     (7 observations deleted due to missingness)     Null Deviance:      159.6      Residual Deviance: 129.4    AIC: 135.4               # LIKELIHOOD RATIO TEST BETWEEN THE 2 MODELS          > J_LRT_TEMP <- lrtest (J_mod.a, J_mod.TEMP)      > J_LRT_TEMP          Model 1: as.factor(Place) ~ as.factor(TempF) + as.factor(AgType) + FemDen     Model 2: as.factor(Place) ~ as.factor(AgType) + FemDen          L.R. Chisq         d.f.            P      2.371378e+01 2.000000e+00 7.089530e-06       I thank you in advance for all of you who can take some of your time to let me know what am I doing wrong, or what am I missing here. Best.


title: visualizing orientation (zenith and azimuth) data (e.g. with a circular or spherical histogram)?

text: I want to understand the distribution of orientation angles in a 3D image (represented as a triangular mesh). For each triangle in the image I have zenith ($\Theta$) and azimuth ($\phi$) ![enter image description here](http://i.stack.imgur.com/APStV.png) And I can compute histograms for each: ![enter image description here](http://i.stack.imgur.com/f5JBg.png) Is there a way that I can combine the information contained in these images into a single histogram representing the joint distribution of these variables?


title: How to generate the time series from a given model?

text: we consider a sparse autoregressive time series of length 1000 obeying the model X(t)=0.2X(t-1)+0.1X(t-3)+0.2X(t-5)+0.3X(t-10)+0.1X(t-15)+Z（t） with nonzero coefficients at lags 1,3,5,10 and 15,where the innovations Z(t) are i.i.d.Garussians with mean zero and standard deviation 0.1. the question is how to simulated 1000 time series from the model with R or SAS? thanks!


title: What to do with collinear variables

text: Disclaimer: This is for a homework project. I'm trying to come up with the best model for diamond prices, depending on several variables and I seem to have a pretty good model so far. However I have run into two variables that are obviously collinear:               >with(diamonds, cor(data.frame(Table, Depth, Carat.Weight)))                        Table       Depth Carat.Weight     Table         1.00000000 -0.41035485   0.05237998     Depth        -0.41035485  1.00000000   0.01779489     Carat.Weight  0.05237998  0.01779489   1.00000000      Table and Depth are dependent on each other, but I still want to include them in my predictive model. I did some research on diamonds and found that Table and Depth are the length across the top and the distance from top to bottom tip of a diamond. Since these prices of diamonds seem to be related to beauty and beauty seems to be related proportions, I was going to include their ratio, say $\frac{Table}{Depth}$, to predict the prices. Is this standard procedure for dealing with collinear variables? If not, what is? Edit: Here is a plot of Depth ~ Table: ![enter image description here](http://i.stack.imgur.com/utIua.png)


title: From relational database to graph: file formatting and software

text: I am approaching graph theory (I have a background in social sciences) to map a discussion forum. My forum/network has four types of vertices: author of thread, thread, comment and author of comments (users writing comments do not coincide with users writing thread). In this case I call "thread" the first post that is then commented. I have scraped the data from the forum website and stored it in a four-table relational database (one table for each type of vertex). What I want to do now is to start playing with the data by visualizing the forum as a graph and compute network properties. But I also need to filter the database for spam comments before analysing the data. The idea would be to describe the forum as different networks (for example by projecting the one set of vertices to obtain a "one-mode" projection). 1) Which file format should I use to dump my database into? I guess that depending on how I built the file I will then have different graphs, but still there is any file format which is particularly handy for manipulating the data and switching from a one-mode to a bipartite graph without recompiling the file from the database? 2) Which software should I contemplate using? (I have some knowledge of R)


title: glm or glmm model with unequal variance

text: I am applying a GLM model with binomial family:               glm(response ~ Treatment, family = binomial, data=dat)      The only explaratory variable _treatment_ is a categorical variable. The model fit turns out that the residuals clearly do not have equal variances among the treatment levels (heteroscedasticy). ![enter image description here](http://i.stack.imgur.com/W8jvT.png) My question is 1) what I can do in this situation? A gam() is not meaningful for categorical variables. On the other hand, gls() is only applicable for linear model. 2) what if I have an extra random factor in the model (GLMM) and the residuals are still heterogeneous?               glmer(response~ Treatment + (1 | groupID), family = binomial, data=dat)      Is there any techniques or R packages can solve this? This question is different from Regression modelling with unequal variance which asks similar situations with a linear model.


title: Intuitive explanation of "integrate out random effect"

text: We are trying to figure out an intuitive reasoning behind integrate out the unobserved random effect. The specific formula is: $f\big(y_i|x_i;\beta, \sigma_c\big)=\int_{-\infty}^{+\infty}\Big(\prod_{i=1}^{T}f_t(y_{it}|x_{it};c,\beta)\Big)(1/\sigma_c)\phi(c/\sigma_c)dc$ I think we get most of the idea: you give a random effect to a person to model the unobserved effect such as IQ or his health. But this random effect is unobserved, so you integrate it out. But here comes my question which makes me doubt my understanding: do you calculate the likelihood for all the possible values the random effect can take? Why do you then say that you have integrated it out? Does it mean that you let the random effect take all possible values and just pick the most probably $\beta$'s, at whatever value of the random effect?


title: How to calculate overlap between empirical probability densities?

text: I'm looking for a method to calculate the area of overlap between two kernel density estimates in R, as a measure of similarity between two samples. To clarify, in the following example, I would need to quantify the area of the purplish overlapping region:               library(ggplot2)     set.seed(1234)     d <- data.frame(variable=c(rep("a", 50), rep("b", 30)), value=c(rnorm(50), runif(30, 0, 3)))     ggplot(d, aes(value, fill=variable)) + geom_density(alpha=.4, color=NA)      ![enter image description here](http://i.stack.imgur.com/gEep0.png) A similar question was discussed here, the difference being that I need to do this for arbitrary empirical data rather than predefined normal distributions. The `overlap` package addresses this question, but apparently only for timestamp data, which doesn't work for me. The Bray-Curtis index (as implemented in `vegan` package's `vegdist(method="bray")` function) also seems relevant but again for somewhat different data. I'm interested in both the theoretical approach and the R functions I might employ to implement it.


title: Creating observed/expected ratio using logistic regression

text: I am using logistic regression to benchmark the performance of some students in different years. I created a scenario as below:               mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")     benchmark.data <- mydata[1:300,] # students form year 1990-1995 as benchmark     compare.data <- mydata[301:400,] # students from year 1996          # logistic regression model created using benchmark student result     temp.glm <- glm(admit~gre+gpa+rank,data=benchmark.data,family="binomial")          # using the regression model to predict how students in 1996 perform     compare.data[,"predict"] <- predict(temp.glm,newdata=compare.data,type="response")          # making a threshold that if the predicted chance of admit > 0.5, then it is asssumed that the student will get admitted     compare.data[,"predict_admit"] <- ifelse(compare.data[,"predict"]>0.5,1,0)     table(compare.data[,c("admit","predict_admit")])          #      predict_admit     # admit  0  1     #     0 59  6     #     1 26  9      From the table, it is seen that 15 students predicted to get admitted and actual number of students get admitted is 35, so the observed/expected ratio is `35/15=2.33`, as it is larger than `1`, so I will say that students in year 1996 is performing better than benchmark. Can I draw my conclusion using the method mentioned above? Besides, how should I set the threshold? Or should I `sum(compare.data[,"predict"])` and treat it as expected value? ### Update 1 I tried and used ROC curve to determine the threshold:               library(ROCR)     benchmark.data[,"predict"] <- predict(temp.glm,newdata=benchmark.data,type="response")     preds <- prediction(benchmark.data[,"predict"],as.numeric(benchmark.data[,"admit"]))     plot(performance(preds,"tpr","fpr"),print.cutoffs.at=seq(0,1,by=0.05))      And the charts suggests that threshold at 0.35 seems to give maximized sensitivity and specificity.


title: Creating obligatory combinations of variables for drawing by random forest

text: **Problem** For my machine learning task, I create a set of predictors. Predictors come in "bundles" - multi-dimensional measurements (3 or 4 - dimensional in my case). The hole "bundle" makes sense only if it has been measured, and taken all together. The problem is, different 'bundles' of predictors can be measured only for small part of the sample, and those parts don't necessary intersect for different 'bundles'. As parts are small, imputing leads to considerable decrease in accuracy(catastrophical to be more accurate) **Possible solutions** I could create dummy variables that would mark whether the measurement has taken place for each variable. The problem is, when random forests draws random variables, it does so individually. So there are two basic ways to solve this problem: 1) Combine each "bundle" into one predictor. That is possible, but it seems information will be lost. 2) Make random forest draw variables not individually, but by obligatory "bundles". **Problem for random forest** As random forest draws variables randomly, it takes features that are useless (or much less useful) without other from their "bundle". I have a feeling that leads to a loss of accuracy. **Example** For example I have variables `a`,`a_measure`, `b`,`b_measure`. The problem is, variables `a_measure` make sense only if variable `a` is present, same for `b`. So I either have to combine `a`and `a_measure` into one variable, or make random forest draw both, in case at least one of them is drawn. **Question** What are the best practice solutions for problems when different sets of predictors are measured for small parts of overall population, and these sets of predictors come in obligatory "bundles"? Thank you!


title: What's more appropriate in the given case: ANOVA or MANOVA?

text: recently, I had an ongoing discussion about ANOVA versus MANOVA, but yet, I am lacking of arguments for MANOVA. My colleque has tracked data for two times of measurement. It's known that there are two groups of people. Both groups are compared regarding their development over the time (t1 to t2) regarding psychometric scaled. At the moment, there are used several ANOVAs, each one per psychometric scale. I am argumenting that a MANOVA fits better because of the fact that multiple ANOVAs lead to cumulation of alpha mistakes. My colleque counter-arguments that the variables are measuring non-dependent constructs wherefore we can justify multiple ANOVAs. Second issue is that there is no option in SPSS offering MANOVA with measurement repetition. What do you think? Best wishes


title: How can I do a chi-square test without being given the critical value table?

text: I have a random number generator that generate integers within `[0, r)`. I want to write a piece of code to test whether the numbers from it are truly uniformly distributed using chi-square test. I select `r` to be very big like `1,000,000`, and I generate random integers `10,000,000` times. Then I get the `X^2`. The next thing I need to do is to check the chi-square value against a degree of freedom table. But I am not allowed to have one, nor using some online distribution calculator. What should I do? Someone told me that I can just calculate an interval: `r^2 - 2*r*(sqrt r)` and `r^2+2*r*(sqrt r)`, then to see whether the value fall in that interval or not. But he doesn't know the reason either. Is the way above correct?


title: Is the F test for R² in (multiple) Regression one- or two-tailed?

text: I have been wondering about the _F_ test that is provided by many statistical packages along with the standard regression output. As I understand it, _F_ can be computed by $$ F_{df_{reg},df_{res}} = \frac{R^2/df_{reg}}{(1-R^2)/df_{res}}. $$ The hypothesis tested by this test can be formulated in two different ways: $H_{0}$: P$^2 = 0$ $H_{1}$: P$^2 > 0$ or $H_{0}$: All $\beta_{i} = 0$ $H_{1}$: One or more $\beta_{i} \neq 0$ The first two hypotheses seem to suggest that the _F_ test is one-tailed, which seems to be inline with my intuition since $R^2$ can not take negative values. The second set of hypotheses, however, suggest two-tailed testing. They, to me, also seem to suggest a direct correspondence between the outcome of the _F_ test for the entire model ($R^2$) and the _t_ tests for the individual coefficients (which I know may not always be the case). So my question comes down to this: **Is the _F_ test, testing whether a model explains a significant amount of variance compared to the null model, a one- tailed or a two-tailed test and why?**


title: Extracting decision function variable from libsvm

text: I'm trying to use LIBSVM's single class SVMs for some classification and need to extract the following sum post classification (i.e. the variable that the decision function takes in) $$ \Sigma_{i=1}^{N} \alpha_i K(z,z_i) $$ I'm not too sure on this so I just want to check that I have the right idea, but could I calculate the sum just by using the primal variable $\textbf{w}$? So for a libsvm model $m$ I just call               w = (m.SVs)'*(m.coef);      and then calculate $\textbf{w}^T \textbf{z}$, where $\textbf{z}$ is the feature vector of what I'd like to classify. I'm not sure though, I guess it may just be easier calculating the kernel function explicitly, but if anyone knows of a good way to access the info in libsvm I would be very grateful


title: Fast Kendall Tau in Matlab (and if necessary another language)

text: As I understand, there are two ways to compute Kendall Tau for a data matrix $X$ (n rows of data points in $\mathbb{R}^p$):   1. corr(X,'type','kendall')   2. http://www.mathworks.com/matlabcentral/fileexchange/27361-weighted-kendall-rank-correlation-matrix/content/kendalltau.m The latter works very, very fast for "small" $X$ (say 300x400). But it runs out of memory when $X$ is say $100\times 1000$. Corr() works very, very slowly in this setting. It can take even an hour sometimes. Are there any better solutions? Ideally in Matlab, but perhaps in another language?


title: Model selection in offline vs. online learning

text: I've been trying to learn more about online learning lately (it's absolutely fascinating!), and one theme that I haven't been able to get a good grasp on is how to think about model selection in offline versus online contexts. Specifically, suppose we train a classifier $S$ offline, based on some fixed data set $D$. We estimate its performance characteristics via cross- validation, say, and we select the best classifier this way. This is what I've been thinking about: how, then, do we go about about applying $S$ to an online setting? Can we assume that the best $S$ found offline will also perform well as an online classifier? Does it make sense to gather some data to train $S$, then take that same classifier $S$ and "operationalize" it in an online setting with he same parameters found on $D$, or might another approach be better? What are the caveats in these cases? What are the key results here? And so forth. Anyway, now that's out there, I guess what I'm looking for is some references or resources that will help me (and hopefully others, who have been thinking about this kinda thing!) make the transition from thinking solely in offline terms, and develop the mental framework to think about the issue of model selection and these questions in a more coherent way as my reading progresses.


title: How good is my error?

text: I'm trying to calculate how good are my measurements in machine learning! Let's say that I have five choices, and that error is 4, 2, 0.002, 3, 6. Naturally, I will pick third one for the hit, but I would like to say following:   * I'm X% certain that hit is third pick    * I'm Y% certain that hit is first (last) pick Of course, X>>Y but I would like to know if is there some kind of statistics metric for this. I can always introduce my metric, but I'm interested is there some already known, so that I can avoid my explanation and just refer to particular metric? The closest one is 'confidence interval', but that gives me the probability that my results are in some interval.


title: Could Bessel's correction make sample variance estimation even more biased?

text: It is well known that Bessel's correction creates an unbiased estimator of variance. What it basically does is divide by $n-1$ instead of $n$. Now what I did is that I chose a few number, like $1,2,3,4,5,60$ and calculated it's population variance which is $452.92$. I then took all possible 4-combinations (15 altogether) and calculated their sample variance (dividing by $n-1$) respectively. The average sample variance is $543.5$ which is off by $90.58$. When I take the population variance (dividing by $n$) of the samples instead I get an average variance of $407.63$ which is off by only $45.29$! I did several other experiments of the same kind with different numbers, population and sample sizes, all with this strange result that the population variance of the samples is less biased than the so called unbiased sample variance. How can that be? What am I missing? **EDIT**   Because of the illuminating discussion in the comments I posted this follow-up question:   Unbiased estimator of variance for samples *without* replacement


title: Rules of thumb for "modern" statistics

text: I like G van Belle's book on Statistical Rules of Thumb, and to a lesser extent Common Errors in Statistics (and How to Avoid Them) from Phillip I Good and James W. Hardin. They address common pitfalls when interpreting results from experimental and observational studies and provide practical recommendations for statistical inference, or exploratory data analysis. But I feel that "modern" guidelines are somewhat lacking, especially with the ever growing use of computational and robust statistics in various fields, or the introduction of techniques from the machine learning community in, e.g. clinical biostatistics or genetic epidemiology. Apart from computational tricks or common pitfalls in data visualization which could be addressed elsewhere, I would like to ask: **What are the top rules of thumb you would recommend for efficient data analysis?** ( _one rule per answer, please_ ) I am thinking of guidelines that you might provide to a colleague, a researcher without strong background in statistical modeling, or a student in intermediary to advanced course. This might pertain to various stages of data analysis, e.g. sampling strategies, feature selection or model building, model comparison, post-estimation, etc.


title: How to plot and analyze a simple grain yield trial?

text: I am looking for a common trial design for comparing 6 different sunflower cultivars, and a way to analyze it as a part of a course I am taking. The common know-how in my class is to set the trial in 4 randomized blocks, each plot with 50-100 plants, gather the grain as bulk from each plot, and compare the means in TUKEY on the jmp software by SAS. Is that the right way to approach this yield comparison? or is there a better way to perform the trial? (A book/article that can guide me step by step to design and analyze this kind of trial will be also helpful) .....


title: Intuition behind which eigenvectors to use in PCA for orthogonal regression

text: I am learning PCA, and while this seems like an obvious question, I can't seem to understand some of the ideas behind which eigenvectors to use when performing orthogonal regression with PCA. My code is below, and a diagram. I am using a parametric model of M=[X.^2;X;Y] to estimate a quadratic with parameters [2;0.7;2]. I am performing eigenvalue decomposition on that. In doing so, I am finding out the directions of variance in the data, and the scale of such variance. PCA aims to find which features best describes the data? So when using the below code, I find that estimating with the first two eigenvectors (directions) (V(:,1:2)') provides a good fit of the function (shown by the green circles in the diagram). However, I cannot get a parameter estimate ([2;0.7;2]) from this data (maybe you can, but I couldn't figure out how). So instead I looked at the final eigenvector (V(:,3)), and this gives me the parameters of my function (by taking V(1:2,3)/V(3,3)). Applying this results in an estimation as shown by the red crosses. As you can see from the diagram, both visual fits are exact. My question is, how can I have estimated this function in two separate ways? I thought that if I used an eigenvector that did not have the smallest eigenvalue, then I would be finding a line in some other direction, not the main direction? Any help anyone can give me would be appreciated, I think I understand how to apply PCA, just unsure as to this directionality thing.               % PCA on Quadratic using parametric model          % Setup data     X=linspace(-1,1,100);     Y=2*X.^2+0.7*X+2;     M=[X.^2;X;Y]';          % Normalize around centroid     meanM=mean(M);     MAdjust = M-repmat(meanM,100,1);          % Perform eigenvalue decomposition     [V D]= eigs(cov(MAdjust));          % Estimate function using V(:,1:2)     FinalData1 = V(:,1:2)'*MAdjust';     MHat1=repmat(meanM,100,1)+(pinv( V(:,1:2))'*FinalData1)';          % Plot starting data, and prediction using V(:,1:2)     figure;hold on;     scatter3(M(:,1),M(:,2),M(:,3),15,'k');     scatter3(MHat1(:,1),MHat1(:,2),MHat1(:,3),50, 'g','d');          % Estimation function using V(:,3)     c = meanM*(V(:,3)/V(3,3));     param = -V(1:2,3)/V(3,3)     yhat =  [X.^2;X;X.^0]'*[param ; c];          % Plot the estimated function     scatter3(X.^2,X,yhat,50, 'r','+');          xlabel('X^2');ylabel('X');zlabel('Y');      ![3D plot of 2*X.^2+0.7*X+2 with predictions from eigenvalue decomposition parameter estimates](http://i.stack.imgur.com/Aymch.jpg)


title: confidence intervals for values estimated from the nonlinear regression model

text: I have a question about nonlinear regression and confidence intervals for values estimated from the model. Here is my problem. I have sets of data where $X$ is the logarithm of the dose of a chemical substance and $Y$ is the response readout (normalized to be between 0 and 100%). I am fitting logistic function to this data using nonlinear regression with 4 parameters: $$y = A + \frac{B-A}{1+e^{-ax+b}}$$ Once I have estimated the $A, B, a,$ and $b$ parameters, I can estimate the value of $X$ at which $Y$ should take 20%, which is what I need. Everything up to this point is clear. However, something is unclear to me: **How do I calculate 95% confidence interval for that estimated $X$ value?** I know there is software out there, which can probably do this for me, but I have to understand the principle and algorithm myself, because I have to implement it in my own library. Any help will be appreciated. **EDIT** Dear Eupraxis and Maarten, thank you for your response: **@Eupraxis** - Here is the picture, the way I understood your explanation about using function prediction intervals to calculate 95% CI bounds of estimated X value. Is it what you proposed (sorry for hand drawings)? ![enter image description here](http://i.stack.imgur.com/cZafb.jpg) What if my function cannot be re-expressed in linear way? I am not even sure how to do it for logistic curve anyway. Please also see my comment for Maarten. **@Maarten** - Thank you for the article! I tried to read it, but frankly it was quite hard to comprehend. I found what it seems to be a much simpler description of delta method for calculation of prediction interval in nonlinear regression: How to compute prediction bands for non-linear regression? Then I was planning to use these 95% PI to calculate 95 CI for x' as Eupraxis suggested above. However, there is one thing I do not quite get in this delta method. Specifically, as described: > First, let's define G|x, which is the gradient of the parameters at a > particular value of X and using all the best-fit values of the parameters. > The result is a vector, with one element per parameter. For each parameter, > it is defined as dY/dP, where Y is the Y value of the curve given the > particular value of X and all the best-fit parameter values, and P is one of > the parameters.) Could you maybe explain it, perhaps with some numeric example I can follow? The rest of the algorithm seems to be clear to me. Thank you for your help in advance, Jul


title: What is PRINCIPAL HESSAIN Direction Model, how and where can I use it?

text: I'm a M.Tech student going through my academic project-work, here I have asked to develop a optimal design and a best equation to fit the data for a Leaching, Solvent extraction and Electrowinning process. So, regarding which I have tried a Uniform Design, simplex Design, Doehlert Design and trying a CCD too. Here the thing is which model I have to use for this kind of design, Regression is not handling with me so how and what I have to do now, please anybody help me. what is Hessian direction and conical Analysis??


title: Fit a GARCH (1,1) - model with covariates in R

text: I have some experiences with time series modelling, in the form of simple ARIMA models and so on. Now I have some data that exhibits volatility clustering, and I would like to try to start with fitting a GARCH (1,1) model on the data. I have a data series and a number of variables I think influence it. So in basic regression terms, it looks like: $$ y_t = \alpha + \beta_1 x_{t1} + \beta_2 x_{t2} + \epsilon_t . $$ But I am at a complete loss at how to implement this into a GARCH (1,1) - model? I've looked at the `rugarch`-package and the `fGarch`-package in `R`, but I haven't been able to do anything meaningful besides the examples one can find on the internet.


title: How to test the mean difference of two Gamma distributions

text: If I don't have the data, but only have the estimates of mean and variance of two independent Gamma distributions. What type of test I can use to test the null hypothesis μ1=μ2?


title: Repeated measures with single measurements

text: I have the following situation: General practitionar (gp), patient (pat) and consultation (cons). Each gp has several patient and each patient can have 1 and more consultations with a specific characteristics cons_x. The outcome is y=0/1 on level consultation. A possible model could be (in lmer/R notation)               y ~ gp_sex + pat_sex + pat_age + cons_x + (1|gp) + (1|gp:pat)      gp is a random variable and pat is nested in gp. Since the outcome is 0/1 (referral to specialist yes/no) the suitable model is a multilevel logistic regression. (1) The problem: 50% of patients have only one consultation, that is, there is no repetition, others have 2, 3 until 6 consultations.I would like to know how to cope with this kind of situation. With other words: I have groups (patients) with only one consultation. So, within the group is perfect correlation due to the only member in the group. Does the mixed models framework cope with this? (2) Another question which arises: The number of consultations weights the effect on patient since the same patient has as many entries as he has consultations. The more consultation a patient has the more patient characteristics are considered. This wouldn't be the case if all patient would have the same number of consultations. Every help and hint is appreciated. Update:   To simplify the analysis I could aggregate on level patient: If a patient has several consultations and at least one consultation has 1 then the patient has outcome 1 otherwise 0. Then I have a simpler model:               y ~ gp_sex + pat_sex + pat_age + (1|gp)      In this case problem (1) will be obsolete but than still I have problem (2) since each gp has different numbers of patients.


title: What methods can be used to transform data?

text: I am solving a binary classification problem with 4 predictor variables. The variables didn't seem to be linearly separable. I have used **Neural Networks and Kernel SVM** which work and give desired accuracy , but in turn are too complex to interpret and have latency issues. Are there any transformation like **box cox or power transform** that I can apply on the data and then use **logistic regression/ decision tree** for classification. I can sacrifice few % of accuracy to get a simple interpretable model. What different methods can be used for data transformation


title: Analyzing changes in a repeated categorical measure

text: I am trying to analyze why some individuals change their response to a categorical survey item Y between time T1 and T2, given some explanatory variables X1 and X2 Specifically, the data looks like this                   ID Y  T  X1 X2            1  1  1  a  b         1  0  2  a  b         2  1  1  b  b          2  1  2  b  b         3  0  1  a  a         3  0  2  a  a              ...      So far I transformed the response variable Y into the four possible combinations of Y at T1 and T2. Specifically,               Y_new  Y_T1 Y_T2      A       0   0      B       0   1      C       1   0      D       1   1      and then I ran a multinomial logit model on the new variable, using the covariates X1 and X2. Is this a feasible approach ? I worry that the transformation violates the IIA assumption. Also, are there any other methods that are better suited for this kind of data?


title: Test for bivariate outliers

text: I have been given a set of data points $(x_i,y_i)$. I have to plot a scatter plot and determine if there are any outliers. But I haven't been taught a method to measure which data point is an outlier and which is not. So how can I do it for example in Sage or R? I found by Google that there is at least two tests to do that, Dixon's and Grubbs's test, so which one should I learn in this problem?               x = c(1,34,6,47,10,49,23,32,12,16,29,49,28,8,57,9,31,10,21,26,31,52,21,8,18,5,18,           26,27,26,32,2,59,58,19,14,16,9,23,28,34,70,69,54,39,9,21,54,26)      y = c(47,76,33,78,62,78,33,64,83,67,61,85,46,53,55,71,59,41,82,56,39,89,31,43,29,55,           81,82,82,85,59,74,80,88,29,58,71,60,86,91,72,89,80,84,54,71,75,84,79)


title: metric on voting data

text: I'm studying computer science but am lost with statistical methods. An algorithm I've been studying needs a formal metric on the space of our data set. Our data set can be described as `voting data'. That is, each data point is a vector containing +1 or -1 as its entries, representing a voter's choices for/against a decision. (We still have to decide how for cases where a voter does not vote, maybe just put a zero on that specific entry.) What kind of distance should I impose on my data for me to obtain a finite metric space on our data set? Using Euclidean distance might be overkill as our vectors do not contain too much geometric information. On the other hand, cosine similarity may not give me a finite metric space. Would the correlation metric (like the one used by _pdist_ in Matlab) suffice?


title: Counting how many times a package has been loaded in R?

text: Following this question, I wish to have some way of counting how many times I am using a package in my daily work. Is there a function/package to do that? In case there isn't, how would you construct such a capability? The way I would do that is by changing it so that at the end of any R session, a log file of the commands would be saved to some location. On that file, I would grep out all instances of "library" and "require". Then save the results into a file (with some time stamp). Lastly, I might want some function to (once in X time) send this file to a remote FTP location - so that other R users could analyse the results. (is there a way to do that with R ?!) If someone wants to try and construct such a machanism - I'd be glad to help by providing a relevant FTP account and by spreading the word on it on "R bloggers" for the good of the community. p.s (mainly for Shane): I wasn't sure if this question should go on stackoverflow or here. If this type of question wasn't debated yet on the meta.stat - it should be. If it was, I'll be glad to know what the conclusion of that discussion was.


title: How build smoothed curve of incidence rates using subject data?

text: I have a dataset with subjects ($\text{data}_0$). Each subject has a follow up time ($\text{fuy}$), an indicator for incidence ($\text{ind} = 0/1$) and a variable which indicates, for example, year of diagnosis ($\text{year}$). One way (I call it ($\text{A}$)) to show the incidence rates according to year would be to compute for each year incidence rate and CI ($\text{data}_1$) and plot them against years, for example, with plotCI (library = gplots).   Another graph would be to plot incidence rate against year and overlay a smoothed curve with CI-region ($\text{B}$). The advantage of ($\text{A}$) is that it considers the standard error for each year whereas this is not the case in ($\text{B}$). ($\text{B}$) considers the noise among all indicidence rates. I wonder if it is possible to plot graph ($\text{B}$) directly from the original subject data ($\text{data}_0$) considering the uncertainty of each year and the jittering during the years for the CI-region (plot ($\text{C}$)). Below is a fake example for clarification. I hope someone can give me some advice. Thanks for help. giordano               # -- create fake data     # n subjects per year     n <- 20     year <- rep(c(1981:2000),each=n)     head(year)     # follow up time in years for each subject     fuy <- 5 + rnorm(400,0,1)     hist(fuy)     sum(fuy<=0)     # indicator of incidence: increasing incidence per year     ind <- c()     for (i in seq(5,14.5,0.5)) {       temp <- i > runif(20,1,20)       ind <- c(ind,temp)     }     ind <- as.numeric(ind)     head(ind)          # build data frame     data0 <- data.frame(year=year, fuy= fuy, ind = ind)     head(data0)     with(data0, table(year,ind, useNA='ifany'))          # -- comupute incidence rates and 95%-CI     # "by hand"-method (poisson regression would give the same incidence rates (ir) and    Wald-CI)     library(sqldf)          data1 <- sqldf(paste("SELECT year, SUM(fuy) fuy, SUM(ind) cases, SUM(ind)/SUM(fuy) ir"                     ," FROM data0"                     ," GROUP BY year"                     )               )     data1          # CI: log(ir) +/- sqrt(1/cases)     data1$cil <- with(data1, exp(log(ir) - qnorm(0.975,0,1)*sqrt(1/cases)))         data1$ciu <- with(data1, exp(log(ir) + qnorm(0.975,0,1)*sqrt(1/cases)))          # -- (A) graph with CI     library(gplots)     with(data1,plotCI(x=year, y=ir, li=cil,ui=ciu))     savePlot(paste(getwd(),"/plot-A.jpeg", sep=""),type="jpeg")          # -- (B) graph with smoothing function     ggplot(data=data1, aes(year,ir)) +       geom_point() +       stat_smooth()     savePlot(paste(getwd(),"/plot-B.jpeg", sep=""),type="jpeg")      ![enter image description here](http://i.stack.imgur.com/L6Syz.jpg) ![enter image description here](http://i.stack.imgur.com/nEteq.jpg)